![]() METHODS FOR USE IN THE PROVISION OF SAND PRODUCTION OF A GEOMECHANICAL RESERVOIR SYSTEM, FOR SAND PR
专利摘要:
methods for use in forecasting sand production from a geomechanical reservoir system, for sand production from a geomechanical reservoir system, and, to operate a geomechanical reservoir system to control sand production from the geomechanical reservoir system. systems and methods are provided for use in forecasting sand production in a geomechanical reservoir system. calculating predictions of sand production may include solving a system of partial differential equations that model the geomechanical reservoir system. systems and methods are also provided for use in operating a geomechanical reservoir system based on the prediction of sand production to control sand production in the geomechanical reservoir system. 公开号:BR112012006055B1 申请号:R112012006055-4 申请日:2010-09-17 公开日:2020-09-15 发明作者:Ricky Howard Dean;Joseph Henry Schmidt 申请人:Chevron U.S.A. Inc; IPC主号:
专利说明:
1. TECHNICAL FIELD [0001] This document refers to systems and methods implemented in a computer for use in the prediction of sand production in a geomechanical reservoir system. This document also refers to systems and methods for use in controlling sand production in a geomechanical reservoir system based on prediction. 2. BACKGROUND [0002] Production in a reservoir system is, in general, the phase that occurs after the development of the reservoir, during which fluids from the reservoir, such as hydrocarbons (oil or gas), are discharged from the reservoir. Sand formation is an occurrence in which solid particles from the formation are produced with fluids from the reservoir. The generic term used to describe small particles from the formation (the rock around a drilling well) that can be produced with the reservoir fluid is "sand." The term "fines" has been used in some literature. The reservoir-forming material generally comprises a type of rock having sufficient porosity and permeability to store and transmit reservoir fluids, such as oil or water. Since sedimentary rocks are porous and are formed in temperature conditions that allow the preservation of fossil remains (from which hydrocarbons are obtained), these are the most common type of reservoir rocks (instead of igneous or metamorphic rocks). Examples of sedimentary rocks include, but are not limited to, conglomerates, sandstones, siltstone, shales, limestone, dolomite, halite (salt), salts, plaster and calcium sulfate anhydride. Sedimentary rocks can include a wide variety of minerals, including, but not limited to, quartz, feldspar, calcite, dolomite and clay group minerals. [0003] Sand production is the migration of formation sand caused by the flow of fluids from the reservoir (such as oil) during production. The production of sand can result from shear rupture or traction rupture within the reservoir formation material. Shear rupture can occur when the pressure in the drilling hole is significantly reduced (such as later in the life of a well), which increases the stress near the drilling well, leading to the formation break. Tensile rupture can occur when the porosity and permeability of the reservoir forming material near the drilling well is significantly deteriorated or when the flow rates are extremely high. Under any condition of tensile rupture, the circulating fluid can exert significant tensile strength forces on individual grains in the formation, which, if excessive, can lead to the breakdown of cementation between individual grains, resulting in rupture by tensile and production of sand. [0004] In many cases, the production of sand may be undesirable, as it can restrict productivity, causes erosion of exploration components, prevents access to the drilling well, interferes with the functioning of downhole equipment and presents significant elimination difficulties. The production of sand increases the diameter of the drilling cavity or the drilling hole, reducing the support around the casing. Consequently, the drilling collapses, the cavity becomes larger and, eventually, the production of the drilling well may cease. If sand production is severe, corrective measures may be necessary to control or prevent sand production entirely, such as filling with artificial gravel or sand consolidation. In extreme cases, massive sand formation can occur, in which the production of sand increases uncontrollably, ultimately causing total erosion of the reservoir material forming the base of the well. [0005] Conventionally, measures were put in place at the beginning of a reservoir development project to avoid, as much as possible, the production of sand in its entirety. For example, unconsolidated formations with multiple production zones can be completed with a coated hole. The coated hole exploration process involves building a barrier system around the reservoir drilling hole to prevent, or delay, the appearance or magnitude of sand formation, including cementing of steel tubes in the drilling well and the use of nets, such as an expandable sand screen technology. These measures, however, can have a negative impact on well productivity. That is, the productivity performance of the well with a coated bore may be much less than that of a well with an open bore. In addition, these measures increase the cost of a reservoir development project, as they require resources, time and manpower for its implementation. [0006] The existence of a certain amount of sand production can help to increase the productivity of a well. Therefore, a method for use in predicting sand production from a geomechanical reservoir system would be useful, as it would allow to determine if, and when, sand production would occur, and to what extent, at the beginning of a project development of a well. The use of these predictions allows the construction and operation of a reservoir in order to produce only a limited amount of sand, substantially sustain a cavity generated in the reservoir and increase the productivity of a well. 3. SUMMARY [0007] As disclosed here, computer-implemented systems and methods to be used in the prediction of sand production from a geomechanical reservoir system are provided. The methods and systems include receiving data indicative of a gradual and massive formation of sand from material within the geomechanical reservoir system; calculate, in a computer system, a value of a critical plastic deformation based on an adaptation of a geomechanical model comprising a hardening model for the data received from massive sand formation; calculate, in a computer system, a value of at least one parameter of a production function based on an adaptation of a geomechanical model, comprising a production function, to the data received from gradual sand formation and use the value of critical plastic deformation and at least one value of one or more hardening parameters; wherein the at least one value of the hardening parameter (s) is calculated based on an adaptation of the hardening model to data indicative of plastic deformation of material within the geomechanical reservoir system; where the hardening model models the behavior of plastic deformation of the material within the geomechanical reservoir system; and where the production function provides for the production of sand from the geomechanical reservoir system. At least one value of the hardening parameter or parameters, the value of the critical plastic deformation or the value of at least one parameter of the production function can be sent to a user interface device, a monitor, a readable storage medium by computer, a local computer, or a computer that is part of a network. [0008] Systems and methods implemented in computer for use in the prediction of sand production from a geomechanical reservoir system can also be provided, comprising: receiving data indicative of plastic deformation, gradual formation and massive formation of sand from material within the geomechanical reservoir system; calculate, in a computer system, a value of one or more hardening parameters based on an adaptation of a hardening model to the received plastic deformation data; where the hardening model models the behavior of plastic deformation of the material within the geomechanical reservoir system; calculate, in a computer system, a value of a critical plastic deformation based on an adaptation of a geomechanical model comprising the hardening model to the data received from massive sand formation; calculate, in a computer system, a value of at least one parameter of a production function based on an adaptation of a geomechanical model, comprising the production function, to the data received from gradual sand formation and use a value of at least one value of the hardening parameters and the value of the critical plastic deformation; where the production function provides for the production of sand from the geomechanical reservoir system. At least one value of the hardening parameter or parameters, the value of the critical plastic strain or the value of at least one parameter of the production function can be sent to a user interface device, a monitor, a storage medium computer readable, a local computer or a computer that is part of a network. [0009] In one aspect of the previous methods and systems, the production function is a function f (x), where f (x) = 0, when x = 0, and f (x) = 1, when x = 1; and where x is a function of the critical plastic deformation. In another aspect of the previous methods and systems, the value of at least one parameter of the production function is a value of an exponent of the production function. The sand production function can be given by the expression: where x = ε / εplim, where ε is a plastic deformation invariant, where gP hm is the critical plastic deformation and where m is the exponent value. [0010] The hardening model can be a modified Bigoni-Piccolroaz model. Plastic deformation data can be obtained from a triaxial test. The gradual formation and massive sand formation can be obtained from a hollow cylinder test. The adaptation of the hardening model to the received plastic deformation data can be achieved by regression. [0011] In a given aspect of the previous methods and systems, the adaptation of the geomechanical model comprising the hardening model to the data received from massive sand formation is obtained by solving a system of partial differential equations that models the geomechanical reservoir system; in which the system of partial differential equations comprises a reservoir flow model and the geomechanical model comprising the hardening model; in which the system of partial differential equations is coupled through a fully expanded Jacobian equation; and in which the resolution of the system of partial differential equations includes solving, simultaneously, in a single time interval, the Jacobian totally expanded based on the data received from massive sand formation. [0012] In another aspect of the previous methods and systems, the adaptation of the geomechanical model comprising the production function to the data received from gradual sand formation is obtained by solving a system of partial differential equations that models the geomechanical reservoir system; in which the system of partial differential equations comprises a reservoir flow model and the geomechanical model comprising the hardening model; in which the system of partial differential equations is coupled through a fully expanded Jacobian equation; and in which the resolution of the system of partial differential equations includes solving, simultaneously, in a single time interval, the Jacobian fully expanded based on the data received from gradual sand formation. [0013] Computer-implemented systems and methods are also provided for use in predicting sand production from a geomechanical reservoir system, comprising: receiving data indicating physical or chemical properties associated with the material within the geomechanical reservoir system; generate predictions of sand production, solving a system of partial differential equations that models the geomechanical reservoir system; in which the system of partial differential equations comprises a reservoir flow model and a geomechanical model of the geomechanical reservoir system; in which the geomechanical model comprises a hardening model; in which a sand production criterion is applied to the geomechanical model; in which the partial differential system is coupled through a fully expanded Jacobian; in which the resolution of the system of partial differential equations includes solving, simultaneously, in a single time interval, the Jacobian fully expanded based on the data received from physical properties; and where the generation is implemented in a computer system. The generated sand production predictions can be sent to a user interface device, a monitor, a computer-readable storage medium, a local computer or a computer that is part of a network. The criterion of sand production can be a critical value of a total deformation, a critical value of a plastic deformation invariant or a maximum effective stress. [0014] In addition, computer-implemented systems and methods are also provided for use in predicting sand production in a geomechanical reservoir system, comprising: receiving data indicating physical or chemical properties associated with the material within the geomechanical reservoir system ; defining a grid comprising a plurality of grid cells; generate predictions of sand production by solving a system of partial differential equations that models the geomechanical reservoir system; in which the system of partial differential equations comprises a flow model of the reservoir and a geomechanical model of the geomechanical reservoir system; in which the geomechanical model comprises a hardening model; in which a sand production criterion is applied to the geomechanical model; in which the system of partial differential equations is coupled through a fully expanded Jacobian equation; in which the resolution of the system of partial differential equations includes solving, simultaneously, in a single time interval, the Jacobian fully expanded based on the data received from physical properties; in which the reservoir model and the geomechanical model are calculated in the grid; and where the generation is implemented in a computer system. The generated sand production predictions can be sent to a user interface device, a monitor, a computer-readable storage medium, a local computer or a computer that is part of a network. The criterion of sand production can be a critical value of a total deformation, a critical value of a plastic deformation invariant or a maximum effective stress. [0015] One aspect of the present disclosure provides a computer system for performing the steps of any of the methods and systems disclosed herein, including the previous systems and methods. The computer system comprises one or more processing units; and one or more memory units connected to one or more processing units, containing one or more memory units containing one or more modules that comprise one or more programs that cause one or more processing units to perform steps including the realization the steps of any of the systems and methods disclosed herein, including the previous systems and methods. In the previous embodiments, one or more memory units may contain one or more modules that comprise one or more programs that cause one or more processing units to perform steps comprising sending to a monitor, an interface device. user, a computer-readable tangible data storage product or a tangible random access memory, a result of the systems and methods. For example, as applicable to the method being executed, the result of the system or method that is sent can be at least a value of one or more hardening parameters, the value of the critical plastic deformation, the value of at least one production function parameter or generated sand production predictions. [0016] Another aspect of the present disclosure provides a computer-readable means of storing a computer program executable by a computer to perform the steps of any of the systems and methods disclosed herein, including the previous systems and methods. A computer program product is provided for use in conjunction with a computer having one or more memory units and one or more processing units, the computer program product comprising a computer-readable storage medium having a program mechanism computer code encoded therein, where the computer program mechanism can be loaded into one or more memory units of the computer and cause the computer processing unit (s) to perform steps comprising performing the steps of any of the systems and methods disclosed herein, including previous systems and methods. In the foregoing embodiments, the computer program mechanism may be loaded into one or more computer memory units of said computer and cause one or more processing units of the computer to perform steps comprising sending to a monitor, a user interface device, a computer-readable tangible data storage product, or a tangible random access memory, a result of the system or method. For example, as applicable to the method being executed, the result of the system or method that is sent can be at least a value of one or more of the hardening parameters, the value of the critical plastic deformation, the value of at least a parameter of the production function or the predictions of sand production generated. [0017] Systems and methods for producing sand from a geomechanical reservoir system are also provided, comprising the operation of the geomechanical reservoir system according to a result of any of the methods and systems disclosed herein, including the systems and methods previous ones. [0018] Systems and methods for operating a geomechanical reservoir system to control sand production from the geomechanical reservoir system are also provided. The systems and methods comprise calculating a value of at least one operational parameter based on an adaptation of a geomechanical model comprising a production function to the indicative data of physical or chemical properties associated with materials within the reservoir system; where the production function provides for sand production from the geomechanical reservoir system; and where the, at least one, value of the operational parameters indicates a condition for the production of stabilized sand from the geomechanical reservoir system; and operate the geomechanical reservoir system according to the value of at least one operational parameter. 4. BRIEF DESCRIPTION OF THE DRAWINGS [0019] Fig. 1 is a flowchart of a method for use in the prediction of sand production from a geomechanical reservoir system. [0020] Fig. 2A represents stress components (ox, σy, Ox) for the daily loading of a sample 200 in a triaxial assay (in this example, Ox = Oz). [0021] Fig. 2B represents the stress component (oy) for a consolidometer test, in which a sample 202 is confined radially to a ring 204 and is only loaded in the y (oy) direction. [0022] Fig. 3 represents an exemplary schematic diagram of a hollow cylinder test. [0023] Fig. 4 shows the results of a hollow cylinder test, plotted as sand production (in cubic centimeters (cm3)) versus time (seconds) and an adaptation to the results using material constants of critical strain value = 0.014 and the exponent (m) = 1. The hollow cylinder sand formation test data were adapted to the two parts of mass formation and gradual sand formation of the curve. [0024] Fig. 5A shows an exemplary schematic diagram of a drilling test. [0025] Fig. 5B shows core samples from a drilling test. [0026] Fig. 6 shows the results of a drilling test, plotted as the production of sand (in cubic centimeters (cm3)) versus time (seconds) and adaptations to the results using different values of the exponent (m). The cylindrical core sample was drilled with a hemispherical end. [0027] Fig. 7 is a flow chart of a method for use in the operation of a geomechanical reservoir system based on the result of a prediction of sand production. [0028] Fig. 8 is a block diagram of an exemplary approach for use in modeling sand production from a geomechanical reservoir system, including a reservoir model and a geomechanical model. [0029] Fig. 9 is a block diagram of an exemplary approach for use in modeling sand production from a geomechanical reservoir system, including a reservoir model, a geomechanical model and a thermal model. [0030] Fig. 10A illustrates an example of a three-dimensional grid for use in calculating models. [0031] Fig. 10B illustrates an example of a two-dimensional grid for use in calculating models. [0032] Fig. 11A shows a graph of the flow surface movement in a Drucker-Prager hardening model with shear hardening. [0033] Fig. 11B shows a graph of the flow surface movement in a Drucker-Prager hardening model with “cap” hardening. [0034] Fig. 12A shows graphs of the flow surfaces in the octahedral (deviatoric) plane for the Drucker Prager model modified to four different values of the deviation (K). [0035] Fig. 12B shows graphs of the flow surfaces in the octahedral (deviatoric) plane for the Bigoni-Piccolroaz model modified for four different combinations of beta (β) and deviation (y) values. [0036] Fig. 13 illustrates an exemplary calculation of sand production using the models. [0037] FIG. 14 illustrates an exemplary computer system for implementing the sand production prediction methods disclosed here. [0038] Fig. 15 shows core samples, including sample dimensions, used for a salt water test (using a hollow cylinder test setup, as illustrated in Fig. 3). [0039] Fig. 16A shows graphs of the triaxial test data (black) and the numerical results (data adaptation) of the simulation (gray) for saltwater saturated core samples. [0040] Fig. 16B shows an enlargement of part of the graph shown in Fig. 16A, where the core sample material becomes almost perfectly plastic. [0041] Fig. 17 shows graphs of the results of the triaxial test (black) and the numerical results of the simulation (gray) for the saturated kerosene samples. [0042] Figure 18 illustrates the boundary conditions applied to the core sample simulations to simulate the load in a hollow cylinder test, where (A) shows the mechanical boundary conditions and (B) shows the boundary conditions of flow applied to the samples in the simulation. [0043] Fig. 19 shows the confinement stress, flow and flow pressure measured and subjected to discretization vs. time for the saturated salt water sample. [0044] Fig. 20 shows a graph of permeability calculated as a function of flow. [0045] Fig. 21 shows a graph of the permeability calculated as a function of the confinement stress. [0046] Fig. 22 shows a graph showing the correspondence between the flows subjected to discretization and modeled as a function of time. [0047] Fig. 23A shows the shape of the flow surface in the deviatoric plane for a Type-Drucker-Prager material model (beta = 0 and deviation = 0). [0048] Fig. 23B shows the shape of the flow surface in the deviatoric plane for a Type-Lade material model (beta = 0 and deviation = 0.95). [0049] FIG. 24 shows a graph of the volume of sand production (cm3) vs. time (seconds) from a measurement of a core sample (black line) and different adaptations of a model to the results using grid sizes of 0.05 mm, 0.01 mm and 0.005 mm. A grid size of 0.01 mm was small enough to capture the details of the model and eliminate grid dependencies. The material constants used were: deviation = 0.5, beta = 0, critical strain value = 0.017 and the exponent = 2. [0050] FIG. 25 shows a graph of the volume of sand production (cm3) vs. time (seconds) from a measurement of a core sample (black line) and adaptation to the results using deviation material constants = 0.5, beta = 0, critical strain value = 0.017 and the exponent = 2. The data were adapted to the two parts of massive formation and gradual sand formation of the curve. [0051] FIG. 26 shows a graph of the volume of water measured and subjected to discretization (liters (L) / min) vs. time (seconds) from the measurement of a core sample (black line) and an adaptation to the results using deviation material constants = 0.5, beta = 0, critical strain value = 0.017 and the exponent = 2. [0052] FIG. 27 shows a graph of effective vs. plastic deformation. hardening (psi) for the salt water test. [0053] FIG. 28 shows a graph of the volume of sand production (cm3) vs. time (seconds) from the measurement of the core samples (black line) and different adaptations to the results in which the deviation and the critical deformation limit are varied, as indicated, and beta and the exponent are kept constant at 0 and 1, respectively. [0054] FIG. 29 shows a graph of the volume of sand production (cm3) vs. time (seconds) from the measurement of core samples (black line) and different adaptations to the results where the deviation = 0.5, beta = 0 and the critical strain value = 0.017, and the exponent value (m) is varied, as indicated. [0055] FIG. 30 shows a graph of the volume of sand production (cm3) vs. time (seconds) from the measurement of core samples (black line) and different adaptations to the results in which the deviation = 0.5, beta = 0 and the exponent = 1, and the limit of critical strain is varied. [0056] FIG. 31 shows horizontal strips, from top to bottom, of a core sample, showing the progression of the cavity. The strips were collected at the end of a sand production test. [0057] FIG. 32 shows numerical simulations of the progression of the cavity over time. Each graph is an axially symmetrical vertical band, in which the gray areas show the location and extent of the cavity and the black areas show intact rock (or spacer). The material constants used in the simulation were: deviation = 0.5, beta = 0, critical plastic strain value = 0.017 and the exponent = 2. [0058] FIG. 33A shows a graph of sand production (cm3) vs. time (seconds), calculated using deviation material constants = 0.5, beta = 0, critical strain value = 0.017 and the exponent = 2. [0059] FIG. 34B shows a numerical simulation of cavity size, calculated for a constant confinement stress of 46 MPa and a constant pressure differential of 0.1643398 MPa from 116675 seconds onwards and using deviation material constants = 0.5, beta = 0, critical strain value = 0.017 and the exponent = 2. The gray area indicates the location of the cavity and the black region indicates intact rock (or spacer). [0060] FIG. 34 shows a numerical simulation of sand production (cm3) vs. time (seconds), calculated for a constant confinement stress of 44 MPa and a constant pressure differential of 1.255517 MPa from 113230 seconds onwards and using deviation material constants = 0.5, beta = 0, strain value critical = 0.017 and the exponent = 2. 5. DETAILED DESCRIPTION [0061] Systems and methods for use in predicting sand production from a geomechanical reservoir system are provided. The systems and methods use measured reservoir properties to generate these sand production predictions, which are useful insofar as they allow, at the beginning of the reservoir development project, to determine if and when sand production can occur and at what rate. [0062] These sand production predictions can be used to determine, at the start of a reservoir development project, the type of exploration techniques that can be implemented to drastically reduce sand production over the life of a well or to allow a certain amount of gradual sand formation (a gradual temporal evolution of sand production), but which does not lead to massive sand formation. For example, decisions can be made regarding the design of the barrier system to be used around the reservoir drilling hole, such as, but are not limited to, sand sieving technology or filling with artificial gravel. Decisions can be made as to whether an unconsolidated formation with multiple production zones should be completed with a lined hole. These predictions of sand production can help to reduce the costs associated with the exploration of the wells, since the decision can be made, depending on the permissible degree of sand production, to install less sand production mitigation equipment than it would be used in the absence of predictions of sand production. Sand production predictions can also be used to make decisions about how to operate the reservoir, such as, but are not limited to, pressure drop, productivity, minimum downhole pressure, production zone temperature and pressure flow of fluid into the drilling well to achieve the desired amount of sand production. In addition, sand production predictions can be used to make a decision about the point in the life of a well where techniques should be used and equipment to mitigate sand production. [0063] The systems and methods of prediction of sand production disclosed here have a physical basis and can be used to predict the phenomena of sand appearance, gradual sand formation and massive sand formation. Sand production prediction systems and methods are deterministic, that is, they are based on physical principles and not just correlations. For example, physically based sand formation criteria are applied to calculations to determine whether sand formation occurs, in what quantity and at what rate. [0064] The operation of a reservoir according to the results of the predictions of sand production can lead to substantial savings, reducing exploration costs. The results of sand production predictions can be used to manage sand formation. The sand formation management can be used, for example, to make decisions regarding the type and cost of sand production prevention exploration equipment to be installed, for example, intubated hole and drilling completions or conventional shored fractures, which produce a limited amount of sand. 5.1 Examples of Sand Production Prediction Systems and Methods [0065] The flowchart of Fig. 1 shows the steps in an exemplary system and method for use in the prediction of sand production from a geomechanical reservoir system. [0066] Step 100. In step 100, data indicating physical or chemical properties of material within the geomechanical reservoir system are received. For example, data indicative of physical or chemical properties of materials can be received in the geomechanical reservoir system during one or more phases of sand formation from material. In the example of Fig. 1, data indicative of gradual sand formation from material in the geomechanical reservoir system are received in step 100A. Gradual sand formation refers to a gradual temporal evolution of sand formation from material in the reservoir. In step 100B, indicative data of massive sand formation from material in the geomechanical reservoir system are received. Massive sand formation refers to a sharp increase in sand production. These data include, but are not limited to, a sand formation rate (in units of sand formation volume over time). Other examples of this data include, but are not limited to, type of formation material, porosity of formation material, permeability of formation material, types of fluid in the reservoir, interstitial pressure, temperature and viscosity of a fluid in the drilling well, production zone temperature, fluid flow pressure in the drilling well, force to withstand fluid flow advance in the drilling well and type of fluid flow in the drilling well. Other examples of this data include, but are not limited to, well depth, oil gravity, oil viscosity, rock type net thickness, current pressure in the reservoir, minimum oil content, oil saturation , the permeability and porosity of the rock, the temperature of the system, the transmissibility of the rock formation, the salinity of the water, the existing fracture system, gas cap, inclination angle, well spacing, receptivity, hydrocarbon composition (HC) , minimum miscibility pressure, pressure ratio, gas saturation, boiling point pressure, critical gas saturation, gas ratio and vertical scan factor. [0067] The data received may be indicative of plastic deformation of material in the geomechanical reservoir system. That is, the data can be one or more parameter values that provide a measure of the plastic deformation of the forming materials. Data indicating elastic behavior of materials can also be received in the reservoir system. Examples of such data are, but are not limited to, Young's modulus, elongation limit, a stress-strain curve for the material, a break resistance limit, strain-hardening behavior, tightening behavior, fracture point . As an example, data indicative of the beginning of the plastic deformation can be received (that is, the point of transition from elastic to plastic behavior). For example, an elongation limit can be used to identify an elastic limit of a material, beyond which additional stress on the material can lead to permanent (plastic) deformation. [0068] The determination of a flow criterion for a material can also be used to determine a beginning of plastic deformation of the material. A flow criterion (which can be presented as a flow surface) can be used to indicate a material elasticity limit (and a beginning of plastic deformation) under different stress combinations. Examples of flow criteria that can be applied to isotropic materials, that is, materials that have uniform properties in all directions, are based on a criterion of a maximum principal stress, a maximum principal strain, a maximum shear stress, an energy total deformation and distortion energy. In the case of a flow criterion based on a maximum principal stress, it can be considered that a flow occurs when the greatest principal stress applied to the material exceeds the uniaxial tensile elongation limit. With a criterion of maximum principal deformation, it can be considered that a flow occurs when the maximum principal stress of the material reaches the deformation corresponding to the flow point during a simple tensile test. In the case of a maximum shear strength yield criterion (Tresca yield criterion), a yield can be considered to occur when the shear strength applied to the material exceeds the shear elongation limit. In a total deformation energy flow criterion, it can be assumed that the stored energy associated with elastic deformation, at the point of flow is independent of the specific stress tensor, so the flow occurs when the deformation energy per unit volume is greater than the strain energy at the elastic limit at simple stress. With the distortion energy flow criterion (Von Mises flow criterion), a flow can be considered to occur when the distortion of the material shape exceeds the flow point for a tensile test. Other examples of flow criteria applied to isotropic materials are the Mohr-Coulomb flow criterion, the Drucker-Prager flow criterion and the Bresler-Pister flow criterion. Examples of flow criteria that can be applied to anisotropic materials, that is, materials whose plastic flow behavior shows directional dependence, include, but are not limited to, Hill's quadratic flow criterion, Hill's generalized flow criterion and the Hosford flow. [0069] Data indicative of elastic behavior, plastic deformation behavior, gradual sand formation or massive sand formation, of a material in the reservoir system can be obtained from tests, such as, but are not limited to, a test triaxial compression test, a triaxial extension test, a uniaxial strain test, a consolidometer test and a hydrostatic compression test. Fig. 2A shows the direction and relationship of the stress components (σx, σY, Oz = σx) that can be applied to a cylinder of material 200 in a triaxial test. Fig. 2B shows the stress component (oY) that can be applied in the y direction to a cylinder of material 202 in a consolidometer test; the material can be confined in a containment ring 204 during compression loading so that no tension is applied in the x, z directions. Either or both tests can be used to provide data indicative of elastic behavior, plastic deformation behavior, gradual sand formation and / or massive sand formation, of a material in the reservoir system. The data can be obtained from tests performed on the material implying different loading paths for the material. The data can also be obtained from a hollow cylinder test, in which a flow test and different combinations of axial, spherical or torsional stress are applied to a hollow cylindrical sample of material. Fig. 3 shows a schematic diagram of a hollow cylinder test setup. Fig. 4 shows data obtained from a hollow cylinder test performed on a material and also shows the starting point of the plastic deformation, the region of gradual sand formation and the region of massive sand formation in the data. The data can also be obtained from a drilling test, in which a flow test is performed on a cylindrical sample of a material, after one end of it has been drilled (that is, subjected to a molded explosive charge). Flow tests can be performed with oil, gas, salt water or any combination of these fluids. Fig. 5A shows a schematic diagram of a drill test configuration and Fig. 5B shows samples of materials that have been subjected to a drill test. Fig. 6 shows data obtained from a hollow cylinder test performed on a material and also shows the starting point of the plastic deformation, the region of gradual sand formation and the region of massive sand formation in the data. [0070] Data can be obtained from tests carried out on one or more core samples taken from a real well site or from a proposed well site. For example, data can be obtained from tests performed on a core sample of material that can be taken from a drilling well. These core samples can be collected in different parts of the reservoir representative of the formation where sand formation can occur. In another example, data can be obtained from tests carried out on a synthetic sample that is created in such a way as to have physical and chemical properties similar to the actual well-forming materials, such as from the well-forming regions where it can be sand formation occurs. [0071] Step 102. Values of one or more hardening parameters are derived based on an adaptation of at least one hardening model to the received data. Hardening models can be used to model the plastic deformation behavior of the material within the geomechanical reservoir system. Hardening models are discussed in Section 5.5.2.4 below. Examples of hardening models include the Drucker-Prager model with shear-hardening or “cap” hardening, the modified Drucker-Prager model with tabular hardening, the modified Bigoni-Piccolroaz model and the Matsuoka-Nakai model. Other examples of hardening models include, but are not limited to, the modified Cam-Clay model, DiMaggio and Sandler's generalized coverage model, the Lade model, the Iwan / Mroz multi-surface model and the plasticity model of continuous surface coverage of Fossum and Fredrich. [0072] Examples of hardening parameters that can be obtained in association with the modified Drucker-Prager model are, but are not limited to, α (a multiplier of a first total voltage invariant), a Drucker-Prager exponent, a constant (F) flow, an effective εpplastic deformation and a deviation (K). Examples of hardening parameters that can be obtained in association with the modified Bigoni-Piccolroaz model are, but are not limited to, the deviation (y), beta (β) and a flow constant (F). [0073] Step 102 can involve several steps, as shown in Fig. 1. In step 102, a hardening model is selected. In a next step, the selected hardening model is adapted to the received data. For example, as illustrated in step 102B, the selected hardening model can be adapted to the received data indicative of the plastic deformation behavior of the material in the reservoir (discussed above in connection with step 100). In step 102C, the values of one or more parameters of adaptation of the hardening model to the data are obtained. [0074] The adaptation of the hardening model to the data received from plastic deformation can be performed using any applicable data adaptation method. For example, adaptation can be performed using a regression method, such as linear regression and non-linear regression. Regression packages that perform data regression adaptation are known in the art. The regression can be performed with limited dependent variables, it can be a Bayesian linear regression, a quantile regression or a non-parametric regression. Adaptation can be performed using a statistical method. Examples of packages that can be used to carry out an adaptation of the hardening model to the received plastic deformation data include, but are not limited to, SAS® (SAS Institute Inc., Cary, NC), MATLAB® (The MathWorks, Inc., Natick, MA), R (accessible via the World Wide Web on the Project R website for Statistical Calculation) and Dap (accessible via the World Wide Web on the GNU Operating System website). [0075] Step 104. In step 104, a value of a critical plastic deformation is calculated based on an adaptation of a geomechanical model comprising a hardening model to the received sand formation data. In one example, sand formation data may be received data indicative of gradual sand formation data or received data indicative of mass sand formation data or a combination of the two. Critical plastic deformation can be a critical value of the plastic deformation of the material that indicates a point at which the material fails, that is, it turns into a granulate without binder to form sand and generates a cavity. For example, critical plastic deformation is a critical value of effective plastic deformation of the material. [0076] An applicable geomechanical model can model the stresses, deformations and / or displacements of isotropic materials, transversely isotropic materials, linear elastic materials, porous materials, solid materials or any combination thereof. As an example, the geomechanical model can model plastic deformation behavior of materials. Examples of geomechanical models are discussed in Section 5.5.2 below (which includes a discussion of hardening models). [0077] The adaptation of the geomechanical model comprising the hardening model for the sand formation data can be performed using any applicable data adaptation method. For example, adaptation can be performed using a regression method, such as linear regression and non-linear regression. In another example, adaptation can be carried out by solving a system of partial differential equations, in which the system of partial differential equations comprises the geomechanical model comprising the hardening model. The system of partial differential equations can also comprise a reservoir flow model (discussed in Section 5.5.1 below). The system of partial differential equations can be coupled through a fully expanded Jacobian, as discussed in Section 5.7 below. The resolution of the system of partial differential equations may include resolving simultaneously, in a single time interval, the fully expanded Jacobian based on the data received (see Section 5.7 below), such as sand formation data. The system of partial differential equations can also comprise a thermal model (discussed in Section 5.5.3 below). A process that can be used to solve the system of partial differential equations is discussed below in Section 5.7. [0078] The adaptation of the geomechanical model comprising the hardening model to the sand formation data can be performed as part of a broader calculation of a system of equations that models the geomechanical reservoir system and which can be used to calculate the stresses , deformations and / or displacements that arise when fluids are injected into or produced from a reservoir, as well as when stresses are applied to the boundaries of a reservoir. In some calculations, a reservoir system model, comprising a geomechanical model (comprising the hardening model), a thermal model and a reservoir flow model, is able to solve systems that include porous flow, heat flow and geomechanics. [0079] In one example, the adaptation of the geomechanical model comprising a hardening model to the data received from massive sand formation can be obtained by solving a system of partial differential equations that models the geomechanical reservoir system, in which the partial differential equations comprises a reservoir flow model and the geomechanical model comprising the hardening model, in which the system of partial differential equations is coupled through a fully expanded Jacobian and in which the resolution of the system of partial differential equations includes solving simultaneously , in a single time interval, the Jacobiana fully expanded based on the data received from massive sand formation. The system of partial differential equations can also comprise a thermal model. [0080] Although step 102 was discussed before step 104, it should be noted that any of these steps could have been carried out in the first place. As illustrated in Fig. 1, any of steps 102 and 104 can be performed at this stage of the flowchart. Steps 102 and 104 can also be performed iteratively, in which a result of executing one step can be used in the other step. For example, if step 102 is performed first, then the value of one or more adaptation parameters in step 102 is used in step 104 to calculate the critical plastic strain. If step 104 is performed first, then the critical plastic deformation calculated in step 104 is used in the adaptation in step 102 to obtain values of one or more hardening parameters. In addition, steps 102 and 104 can be performed repeatedly and iteratively to achieve better data adaptations. As a non-limiting example, if the hardening model is a modified Bigoni-Piccolroaz model, then one can obtain an initial value of one or more of the deviation (y), beta (β), α or the constant (l “ ) flow in step 100 is used in step 104 to calculate, for example, a value of a critical plastic deformation that adapts to the indicative data of massive sand formation. The calculated value of the critical plastic deformation can be used in a second calculation iteration in step 100 of a hardening model parameter. The second iteration of step 100 could provide a better adaptation to the data received. The result of the second iteration of step 100 can be provided to step 104 for a second iteration of calculating the critical plastic deformation to improve adaptation to the massive sand formation data. Iterations may cease after the adaptations in steps 100 and 104 converge to better match the data. [0081] Step 106. A value of at least one parameter of a production function is calculated based on an adaptation of a geomechanical model comprising the production function for the data received from gradual sand formation and using the deformation value critical plastic (discussed in step 104) and a value of at least one of the hardening parameters (discussed in step 102). [0082] The production function models the amount of sand production from a material before reaching critical plastic deformation. Production functions are discussed in Section 5.5.2.5 below. The production function can be any function f (x) whose values vary from f (x) = 0, when x = 0, to f (x) = 1, when x = 1. The sand production function can also be any function f (x) whose values can be scaled to be in an interval of f (x) = 0, when x = 0 and f (x) = 1, when x = 1. In another example, the sand production function it can also be any function f (x) whose values can be transformed to fall within the range of f (x) = 0, when x = 0, and f (x) = 1, when x = 1, by applying a suitable transform, such as, but not limited to, an undulate transform and a Lagrange transform. The term x can be any function g () of the critical plastic deformation (that is, x = εp g (εplim)). The function f (x) can be any monotonic function of x, including, but not limited to, a fraction function, a power function, a sine function, a cosine function, a logarithmic function, an exponential function, a sigmoid function or any combination thereof. In some examples, x may be a function of the critical plastic deformation (εplim) and the plastic deformation invariant (εp) of a material, such as, but not limited to, a ratio of x = εp / εplim. In one example, the production function can be given by f (x) = (εp / εplim) m] where m is an exponent. [0083] The parameter for which a value is calculated in step 106 can be any parameter that can be used to characterize the production function. For example, if the production function is a power function, the parameter can be at least an exponent of the power function. In a given P εp example where the production function is given by f (x) = (εp / εplim) m, the parameter can be the exponent m. In the previous example, the at least one parameter of the production function for which a value is calculated in step 106 is an exponent. In other examples, the parameter can be a multiplier of a production function term. [0084] The production function can be used to predict the amount of sand production from a material before failure, that is, during the gradual formation of sand before reaching critical plastic deformation. [0085] For example, Fig. 4 shows the results of an adaptation of the production function p c-P (εp / εplim) m to the data received from a hollow cylinder test. The data were adapted to parts of massive sand formation and gradual sand formation of the curve using critical material deformation value constants εplim) = 0.014 and the exponent (m) = 1. In another example, Fig. 6 shows the p pP result of different adaptations of the production function (εp / εplim) m to the data received from a drilling test using different values of the exponent (m). [0086] The geomechanical model comprising the production function can be adapted to the sand formation data using any applicable data adaptation method. For example, adaptation can be performed using a regression method, such as linear regression and non-linear regression. In another example, adaptation can be performed by solving a system of partial differential equations, in which the system of partial differential equations comprises the geomechanical model comprising the production function. The system of partial differential equations can also comprise a reservoir flow model. The system of partial differential equations can be coupled through a fully expanded Jacobian and the resolution of the system of partial differential equations can include solving simultaneously, in a single time interval, the fully expanded Jacobian based on the received data, such as the data sand formation. The system of partial differential equations can also comprise a thermal model. [0087] In one example, the adaptation of the geomechanical model comprising the production function to the data received from gradual sand formation can be obtained by solving a system of partial differential equations that models the geomechanical reservoir system, in which the system of partial differential equations comprises a reservoir flow model and the geomechanical model comprising the production function, in which the system of partial differential equations is coupled through a fully expanded Jacobian and in which the resolution of the system of partial differential equations includes solving simultaneously, in a single time interval, the Jacobiana fully expanded based on the data received from gradual sand formation. [0088] The adaptation of the geomechanical model comprising the production function to the sand formation data can be performed as part of a broader calculation of a system of equations that models the geomechanical reservoir system and which can be used to calculate the stresses , deformations and / or displacements that arise when fluids are injected or produced from a reservoir, as well as when stresses are applied to the boundaries of a reservoir. In some calculations, a reservoir system model, comprising a geomechanical model (comprising the production function), a thermal model and a reservoir flow model, is able to solve systems that include porous flow, heat flow and geomechanics. [0089] In step 108, information can be sent that can be used to predict sand production from the geomechanical reservoir system. This information can be, but is not limited to, at least one value of the hardening parameter (s), the value of the critical plastic deformation and / or the value of at least one parameter of the production function. Information can be sent to a user or to various components, such as a user interface device, a computer-readable storage medium, a monitor, a user-accessible local computer or a user-accessible computer that is part of a network. For example, the output can be visually presented to a user using a monitor or user interface device, such as a portable graphical user interface (GUI), including a personal digital assistant (PDA). 5.1.1 Other Sand Production Prediction Systems and Methods [0090] Other exemplary systems and methods for use in predicting sand production from a geomechanical reservoir system include a step of applying a sand production criterion to a geomechanical model that comprises one or more hardening models. The system and method comprise the steps of receiving data indicative of physical or chemical properties associated with materials within the geomechanical reservoir system and generation of sand production predictions by solving a system of partial differential equations that models the reservoir system geomechanical. In addition to the geomechanical model, the system of partial differential equations can comprise a reservoir flow model and / or a thermal model of the geomechanical reservoir system. [0091] One or more criteria of sand production can be applied to the geomechanical model. The criteria for sand production can be determined as when a critical value is reached for: (1) the total deformation invariant or (2) the plastic deformation invariant or (3) the maximum effective stress. When the critical value of the sand production criterion is reached, the material fails, that is, it transforms into a granulate without binder in order to form sand and generate a cavity. The criteria for sand production are discussed in Section 5.5.2.5 below. [0092] In one example, the system of partial differential equations can be coupled through a fully expanded Jacobian, in which the resolution of the system of partial differential equations, such as using a computer system, includes solving simultaneously, in a single interval time, Jacobiana fully expanded based on the data received. The generated sand production predictions can be sent to a user, a user interface device, a monitor, a computer-readable storage medium, a local computer or a computer that is part of a network. For example, generated sand production predictions can be presented visually to a user using a monitor or a user interface device, such as a portable graphical user interface (GUI), including a personal digital assistant (PDA). 5.2 Systems and Methods for Operating a Geomechanical Reservoir System [0093] Systems and methods for use in the control of sand production from a geomechanical reservoir system during operation are also disclosed. The method operates the geomechanical reservoir system according to a result of the implementation of any of the sand production prediction systems and methods disclosed here. The flowchart in Fig. 7 shows steps in an exemplary system and method for use in operating a geomechanical reservoir system based on the result of a prediction of sand production. [0094] In step 700, data indicative of physical or chemical properties associated with materials within the reservoir system are received. The data received in step 700 can include any of the data described in step 100 above. [0095] In step 702, a value of at least one operational parameter is calculated based on an adaptation of a geomechanical model comprising a production function to the received data, in which the production function provides for sand production from the system geomechanical reservoir and where the at least one value of the operational parameter indicates a condition for the production of sand from the geomechanical reservoir system. As an example, the calculated value of the operational parameter indicates a condition for the production of stabilized sand from the geomechanical reservoir system. [0096] The prediction of sand production from a reservoir of any of the systems or methods disclosed herein can be used to calculate or obtain a value of at least one operational parameter for operating a reservoir to achieve the desired quantity of sand production or the desired behavior of sand production. Examples of operating parameters include, but are not limited to, pressure drop, productivity, minimum downhole pressure, production zone temperature and fluid flow pressure in the drilling well, confinement stress and pressure differential. The results of the implementation of sand production predictions can be used to calculate the value of one or more operational parameters that lead to the development and maintenance of a substantially stabilized cavity in the reservoir. A cavity in a reservoir can be substantially stabilized if, over time, the cavity does not grow substantially or grow insignificantly. For example, values of one or more operational parameters can be obtained based on results of production sand predictions that indicate an operational condition for the reservoir through which a limited amount of sand is produced from the reservoir and a cavity in the reservoir created by this limited amount of sand formation is substantially stabilized. Sand production from a reservoir can be kept to a minimum by controlling production variables, such as, but are not limited to, pressure drop, minimum downhole pressure and productivity. [0097] In step 704, the geomechanical reservoir system is operated according to the value of, at least one, operational parameter. That is, the geomechanical reservoir system can be operated with the calculated value for at least one operational parameter, such as a pressure drop value, productivity, minimum downhole pressure, production zone temperature, flow pressure of fluid in the drilling well, confinement tension, pressure differential or any combination of these parameters. [0098] The results of implementing any of the sand production prediction systems and methods can be used to determine the type of exploration techniques that can be installed in a reservoir to achieve the desired amount of sand production throughout the life of a well. For example, the design of the barrier system to be used around the reservoir drilling hole (such as, but not limited to, sand screen technology or artificial gravel filling) can be selected based on these results. The predictions of sand production can also influence the exploration strategy of the well, such as open hole, coated hole and drilled coated hole, use of sieve or “frac-and-pack” coatings (hydraulic fracturing followed by injection of an agent fracture). In addition, sand production predictions can be used to make a decision about the point in the life of a well where techniques should be used and equipment to mitigate sand production. 5.3 Examples of Modeling Methods [0099] Figs. 8 and 9 represent examples of systems for use in predicting sand production from a geomechanical reservoir system during oil production. An applicable modeling system includes one or more models that describe various physical aspects of a geomechanical reservoir system. Fig. 8 and Fig. 9 represent modeling systems that include a reservoir fluid flow model and a geomechanical model. The reservoir fluid flow model describes, for example, porous flow, production and injection. The geomechanical reservoir model describes, for example, stresses, strains and displacement that arise when fluids are injected into or produced from a reservoir and when stresses are applied to the boundaries of a reservoir. The system in Fig. 9 also includes a thermal model. The thermal model described a heat flow. A system of non-linear partial differential equations can be used to interrelate the various aspects of these models. [00100] After receiving data indicative of physical or chemical properties of material within the geomechanical reservoir system (such as, but not limited to, plastic deformation, gradual sand formation or massive sand formation), a resolving agent generates predictions (for example, sand production predictions) by applying the steps of the method in Section 5.1 above and, at relevant points, solving the system of partial differential equations. In the resolving agent of Figs. 8 and 9, the system of partial differential equations can be coupled through a fully expanded Jacobian. The resolution of the system of partial differential equations includes solving, simultaneously, in a single time interval, the Jacobian fully expanded based on the data received. [00101] As shown in Figs. 8 and 9, the sand production prediction steps (steps 100 to 106 discussed above) can be performed iteratively with the solution of the system of non-linear partial differential equations. That is, in the stages of prediction of sand production that include the realization of an adaptation, for example, each of the steps 102, 104 and 106, the calculation can be carried out simultaneously solving, in a single time interval, a Jacobian totally expanded form of the equations representing the geomechanical reservoir system that are relevant to the particular stage. The generated sand production predictions can be sent to various components, such as sent to a user interface device, a computer-readable storage medium, a monitor, a user-accessible local computer or a user-accessible computer making part of a network. [00102] The system of non-linear partial differential equations comprises equations that correspond to the models that will be used in the analysis of the geomechanical reservoir system in the specific stage of sand production prediction. For example, Fig. 8 provides an example in which the system of non-linear partial differential equations includes equations corresponding to a reservoir flow model and a geomechanical model of the geomechanical reservoir system. Depending on the sand production prediction step in progress, the geomechanical model may include one or more hardening models, a sand production model (comprising the production function) or both. In another illustration, Fig. 9 provides an example in which the system of non-linear partial differential equations includes equations corresponding to a reservoir flow model, a geomechanical model and a thermal model of the geomechanical reservoir system. Examples of equations that correspond to each of the different models of the geomechanical reservoir system are provided in Section 5.5. [00103] The coupling of the various aspects of the models can be implemented as through variables in a fully expanded Jacobian. For example, a fully expanded Jacobian can function as a coupling of fluid flow in the reservoir to the geomechanical model by one or more of the following variables: effective stress, porosity and one or more displacements associated with the geomechanical model. A fully expanded Jacobian variable that couples the geomechanical model to fluid flow can be porosity and permeability that is associated with the reservoir flow model. A fully expanded Jacobian variable that couples the thermal model to the geomechanical model can be a thermal stress associated with the thermal model. A fully expanded variable in Jacobiana that couples the thermal model to the flow model of the reservoir can be fluid viscosity, conduction and convection in the reservoir associated with the thermal model. The fully expanded Jacobian can include terms related to a rate of change (ie, a derivative with respect to time), a spatial derivative, or partial spatial derivative, of a coupling variable, where the derivatives can be of any order, for example example, a first order derivative, a second order derivative, a third order derivative, etc. First, second, third and / or higher order derivatives (either time-related or spatial derivatives) of the coupling variables can be included in fully expanded systems of equations. Examples of variables that can couple the different models are provided in Section 5.6. [00104] The fully expanded Jacobian system of non-linear equations can be solved using numerical approaches, such as the approach discussed in more detail in Section 5.7, in which the system of non-linear equations is solved, for example, using an expansion Newton-Raphson method completes all solution variables, which increases the solution's stability and allows second-order convergence rates for non-linear iterations. Examples of devices and software implementations of the different methods disclosed here are discussed in Section 5.8. [00105] In another aspect, a system and method may include the steps of receiving data indicative of physical or chemical properties associated with the material within the geomechanical reservoir system, defining a grid comprising a plurality of grid cells and generating predictions of sand production by solving a system of partial differential equations that models the geomechanical reservoir system. In addition to the geomechanical model, the system of partial differential equations can comprise a reservoir flow model and / or a thermal model of the geomechanical reservoir system. [00106] One or more criteria of sand production can be applied to the geomechanical model. The criteria for sand production can be determined as when a critical value is reached for: (1) the total deformation invariant or (2) the plastic deformation invariant, or (3) the maximum effective stress. The criteria for sand production are discussed in Section 5.5.2.5. [00107] In calculations using criterion (2), that is, calculations involving the calculation of the critical value of the plastic deformation invariant (the critical plastic deformation), steps 100-106 can be implemented and it can be calculated at least one parameter of a production function. The production function can be used to predict the amount of sand production from a grid cell before reaching the critical plastic strain value. [00108] In one example, the system of partial differential equations can be coupled through a fully expanded Jacobian, in which the resolution of the system of partial differential equations, such as using a computer system, includes solving simultaneously, in a single interval time, Jacobiana fully expanded based on the data received. The reservoir model, thermal model and geomechanical model can be calculated in three-dimensional grid cells or two-dimensional grid cells. Three-dimensional and two-dimensional grid cells that can be used in the simulation methods here are described in Section 5.4. [00109] The generated sand production predictions can be sent to a user, a user interface device, a monitor, a computer-readable storage medium, a local computer or a computer that is part of a network. [00110] Solving the system of non-linear equations, implicitly, for example, using a complete Newton-Raphson expansion of solution variables, can improve numerical stability (for example, when dealing with cavity generation or any simulation involving very small grid blocks). The use of a Jacobian fully expanded from an implicitly coupled system of equations provides more stability for the solution process. 5.4 Simulation Method [00111] Fig. 10A illustrates an example of a three-dimensional (3D) grid that can be used for calculations of the geomechanical model and the reservoir model and / or thermal model. For example, one or more of the multipoint flow model (such as for the reservoir and porous flow), geomechanical model (such as for the calculation of stress, displacement and / or the generation of the cavity (transformation into granules without binder)) and thermal model can be calculated on the 3D grid. The calculation of the geomechanical model and the reservoir model and / or thermal model can use a parallel processing approach to couple the grids. The dimensions of the grid cells can be in the order of feet, inches or fractions of an inch. In another example, the dimensions of the grid cells can be in the order of meters, centimeters, millimeters, microns or micron fractions. [00112] The 3D grid can be structured or unstructured hexahedral grids comprising hexahedral elements. The hexahedral grid cell has eight corners, twelve edges (or sides) and six faces. Each cell in the hexahedral grid can include at least eight nodes (one in each corner) or more and up to twenty-seven (27) nodes (that is, a node in the center of each face, in the center of each side, in the center each edge and in the center of the cell). Different hexahedral cells can include a different number of nodes. For example, the 3D grid can include structured or unstructured tetrahedral grid elements. In another example, the 3D grid may include other element morphologies that are in the interval between the tetrahedral and hexahedral grid elements, structured or unstructured. The 3D grid can include any combination of the grid elements mentioned above. [00113] A two-dimensional (2D) grid can also be used for calculations of the geomechanical model and the reservoir model and / or thermal model. For example, a 2D grid can be used for axisymmetric calculations. Fig. 10B illustrates an example of a two-dimensional (2D) grid that can be used for calculations of the reservoir model, thermal model and geomechanical model. [00114] The 2D grid can be structured or unstructured quadrilateral grids comprising quadrilateral elements. Each quadrilateral grid cell has four corners and four edges. Each quadrilateral grid cell includes at least four nodes (one in each corner) and up to five nodes (that is, one node in the center). [00115] In certain examples, calculations can be performed on a 2D grid and a 3D grid. For example, depending on the topography of the reservoir, some of the calculations can be performed on a 2D grid, while other calculations can be performed on a 3D grid. In these calculations, the nodes of the 2-D network can be configured to coincide with the nodes of one of the outer borders of the 3D grid and certain calculations, such as fracture widths and 3D displacements, can be coupled at common node points. The input data format for the 2D network can be similar to the format for the 3D grid. [00116] For certain calculations, parameters, such as fluid flow, displacements, cavity generation and traction, can be monitored through the elements of the grid cells. 5.5 Models [00117] Examples of the differential equations that correspond to each of the different models of the geomechanical reservoir system are provided below. The differential equations for the models included in a calculation can be combined to produce a fully coupled implicit formulation. A consistent set of units can be used for all variables in the equations included in a calculation. 5.5.1 Reservoir Model [00118] The system of equations for porous flow includes conservation of mass where <p is porosity and p is fluid density which can be a function of pressure. The model allows exploration of wells in the reservoir elements and qw in the above equation represents the injection in the reservoir elements. [00119] The speed at Darcy's speed relative to the porous material and can be defined by where K is a tensor permeability, p is the viscosity that can be a function of pressure, p is fluid pressure and pgVh is a gravitational term. [00120] The geomechanical variables included in the fluid flow equations enhance the coupling between the flow and strain models (the definitions of certain geomechanical terms are described in Section 5.5.2 below). [00121] Temperature-dependent properties of water can be introduced in several ways in the case of calculations involving changes in temperature. Water properties can be introduced as pressure (P) and temperature (T) functions for fully coupled calculations. In the case of iteratively coupled calculations, the properties of the water can be introduced as pressure functions and then, modification factors for the effects of temperature can be used. The treatment for temperature dependent properties is explained in more detail in Section 5.5.3. [00122] The thermal behavior of fluids can also be modeled by modifying fluid properties using modification factors (described in Section 5.5.3.3 below). 5.5.1.1 Multiphase Porous Flow [00123] The reservoir model allows several models of phase behavior, ranging from single-phase compositions to black oil compositions based on (phase) fugacity. The Darcy flow can be modeled for aqueous phases, non-aqueous liquid phases and non-aqueous vapor phases and nc components. Any models of phase behavior can be used with the porous flow models disclosed here. Fluid flow equations can be presented in terms of a general compositional formulation. Partial differential equations representing component mass balances for multiphase flow are: where Nic is the concentration of ic component per unit of pore volume, given by is the molar fraction of the ic component in phase a, pa is the molar density of the phase, and qiC is the molar flow of the ic component per unit volume of the reservoir. The phase velocity α is given by [00124] Phase pressures can be defined by where Pca is the capillary pressure and P is the reference pressure. The reference pressure is used for PVT calculations, well calculations and geomechanical calculations. The reference pressure can be a non-aqueous phase pressure for two-phase models and the non-aqueous liquid phase pressure for three-phase models. [00125] Porosity can be defined as for porous flow models, where cpo, Cr and the initial pressure Po are location functions. [00126] In the case of calculations that include the geomechanical model, the porosity relative to the initial undeformed apparent volume is given by Equation 15, where it is seen that Cr may be related to the Biot constant for porosity. 5.5.1.2 Non-Darcy Law Flow [00127] In certain examples, the flow not conforming to Darcy's law can be modeled using the Forschheimer equation to modify the relationship between the pressure gradient and the fluid velocity. In another example, flow that does not conform to Darcy's law can be modeled by specifying a general relationship between fluid velocity and pressure gradient. [00128] The Forschheimer equation for velocity not in conformity with Darcy's law (which would replace Equation 2 above), for systems involving 3-D multiphase flow, in anisotropic media, can be given by: where pa is the viscosity of phase a, K is the permeability tensor, kra is the relative permeability of phase a, va is the Darcy velocity of phase α, the expression is the norm defined by and the parameter βk is a coefficient not according to Darcy's law. The βk parameter can vary across the reservoir, therefore, a value for βk can be entered for each grid block. The βk coefficient not according to Darcy's law is related to the inverse of the transition constant, that is, βk = 1 / F. See Barree et al., “Beyond Beta Factors, A Complete Model for Darcy, Forchheimer, and Trans-Forschheimer Flow in Porous Media,” SPE 89325, annual SPE Technical Conference and Exhibition, Houston, TX 26-29 September, 2004. The Reynolds number for the phase is given by the equation [00129] The units of the terms in equation 8 must be chosen so that the result is dimensionless. After combining Equations 7 and 8, the Forschheimer equation becomes in which the non-Darcy law flow is expressed as a phase-dependent permeability modification function that varies with the Reynolds number for that phase. This is different from the standard permeability modification function, because this formulation that does not conform to Darcy's law has separate values for each phase. [00130] In some cases, the Forschheimer equation does not provide an adequate approximation for a flow that does not conform to Darcy's law. Modification functions of the following formula can be used for outflow that does not conform to approximate Darcy's law: [00131] Modification functions can be built that satisfy the restrictions [00132] In the case of the standard Forschheimer equation, the following function can be specified: [00133] For an extended Forschheimer equation, the following functions can be specified: [00134] Equation 3 in J.L. Miskimins JL, et al. “Non-Darcy Flow in Hydraulic Fractures: Does It Really Matter ” article SPE 96389, SPE's annual Technical Conference and Exhibition, Dallas, TX, October 9-12, 2005, is another form of a formulation not in accordance with Darcy's law that is applicable to the methods disclosed in this application. 5.5.1.3 Calculation in the Reservoir Model [00135] The calculation of the reservoir model can be in the 3D grid (which can be the grid used for the geomechanical model). The 3D reservoir grid can include porous flow calculations. Fluid velocity terms can be calculated for flow between cells in the reservoir, as well as flow between cells in the reservoir and cells in the cavity (cells where sand has formed). A major variable for the porous flow equations can be fluid pressure or fluid composition, which can be evaluated at the center of each hexahedral (cell-based) element. In certain calculations, a multipoint flow algorithm can be used for unstructured reservoir flow calculations so that the resulting computational model can be 27 points for general hexahedral elements of a 3D grid when eight elements share a common corner. 5.5.2 Geomechanical model 5.5.2.1 Poroelastic materials [00136] A linear relationship (small strain) can be used for strain-displacement relationships. The coupled flow / displacement model relating stress, strain, temperature and interstitial pressures can be based on Biot's poroelastic theory. The equilibrium equation can be based on total stresses and assumes a quasi-static equilibrium. [00137] Poroelastic equations can be formulated in terms of total stresses, apparent strains, temperatures and interstitial pressures. The total stress can be defined by the average stresses that could be observed in a flat section of the reservoir, where the flat section includes loads carried by the fluid and interstitial pressures of the fluid. Apparent deformations can be deformations that could be observed from a deformation meter if it were connected to the porous material in deformation. [00138] The system of equations for linear poroelastic displacements includes the strain-displacement equation where commas imply differentiation, ui is the displacement in the i direction, εy is the apparent deformation of the porous material and the expansion corresponds to positive deformations. Total stresses satisfy equilibrium equations where the tension tensor is symmetric and the term fi of gravity is a function of solids density, fluid densities and porosity. Traction or displacement boundary conditions can be specified in all three directions in all six boundaries of the three-dimensional grid on which the model is calculated. [00139] When temperature differences are not taken into account, the constitutive equations relating total stresses, apparent strains and interstitial pressure are where the stress is positive, the repeated kk index implies summation, is the initial in-situ stress, p0 is the initial pressure, E is the modulus of elasticity, v is the Poisson's ratio, α is the Biot constant in equations of stress / strain, õij and 1 when i = je 0 when i ψj. Deformations can be assumed to be null when σ ij = σ ij. [00140] In the examples in which temperature differences are taken into account, the constitutive equations are: where CIT is the coefficient of thermal volumetric expansion for stress / strain equations and K is the elastic compressibility module. The pressure po is the initial interstitial pressure and To is the initial temperature. [00141] If the terms of stress and pressure are combined to form σij = σij + αPδij, then the equation becomes a standard, thermal, linear and elastic constitutive equation in which the stresses have been replaced by the effective stresses. If the initial in situ stresses and initial interstitial pressure are zero, then the equation takes the standard form [00142] The relationship between porosity (in relation to the apparent undeformed volume), and the deformations and fluid pressure (when temperature differences are not taken into account) is given by: where Equation 15 assumes that the initial deformations are zero, q> o is the initial porosity and M '1 is the Biot constant for interstitial pressure in the porosity equations. [00143] When differences in temperature and deposition fraction are taken into account, the porosity relative to the apparent undeformed volume is defined as: where α and M 1 are constant Biot, P is the phase pressure (for multiphase flow), v is the coefficient of thermal volumetric expansion for porosity and o is the deposition fraction (the volume fraction of solid waste deposited by volume apparent value, for example, of a grid element in a calculation). Solid waste can be deposited within pores when waste moves through the reservoir and there can be reductions in porosity and permeability when waste deposits accumulate in interstitial spaces. The porosity in a calculation grid element can be a function of fluid pressure, temperature and deformations, while defining that the amount of porosity reduction due to the deposition of solids can be equal [00144] For an isotropic material, the six parameters of poroelastic materials: E, v, α and M'1, OT and av, are determined before the application of geomechanical equations to the modeling of a geomechanical reservoir system. [00145] In certain examples, reservoir permeabilities can be expressed as an initial directional permeability (Kabs) multiplied by a permeability multiplier f for permeability at each point and for each time step: K = Kabsxf, where f is a function of one or more other parameters, such as fluid pressure, total stresses, apparent volumetric deformation, interstitial pressure, initial reference pressure, main stresses, effective plastic deformation, current porosity, initial porosity and deposition fraction. [00146] In the constitutive equations for a transversely isotropic material, deformations can be expressed in terms of stresses. In certain calculations, effective stresses can be used for the calculations and the initial in-situ stresses may be different from zero. In other calculations, the equations can be simplified using total stresses and using an assumption that the initial in-situ stresses can be zero. The constitutive equations relating stresses and strains for a transversely isotropic material can be expressed as: [00147] Equations 17-22 assume that the axis of symmetry is the z direction and that the vertical direction is the z direction. Five elastic constants in Equations 17-22 can be supplied to the model before performing a calculation involving transversely isotropic solids. It can be assumed that the Biot constants are isotropic and that the thermal expansion coefficients are isotropic when performing poroelastic calculations, so that two Biot constants and two thermal expansion coefficients can be provided in addition to the five transversely isotropic constants when analyzing a transversely isotropic porous material. In the case of a poroelastic calculation including the thermal model with transversely isotropic elastic constants, the stresses in Equations 17-22 can be effective stresses and the stresses in Equations 17-22 can be effective stresses defined as: where εy is the true strain. 5.5.2.2 Poroplastic Materials [00148] A poroplastic material exhibits a non-linear behavior, in that it can undergo permanent volumetric deformations (that is, plastic) and, therefore, the porosity changes. Large fluid pressures can lead to the flow of poroplastic material. Consequently, geomechanical calculations for a poroplastic material can predict large and sudden changes in fluid porosities. These large and sudden changes in porosity can cause significant stability problems and also produce negative fluid pressures. Negative pressures usually arise when the compressibility of the fluid is low, the permeability is low and the expansion of porosity is sudden and large. The equations for a poroelastic material discussed in Section 5.5.2.1 above are also applicable to poroplastic materials. However, porosity can be modified to adopt changes in porosity that can be predicted for poroplastic materials. [00149] In certain examples, an equation can be used to cushion sudden changes in porosity, in order to improve the numerical stability of the calculations and to reduce the frequency of encounters with negative pressures. In these calculations, the porosity in the reservoir model can be defined as a fluid porosity (ψfiuid) and can be treated differently from the porosity in the geomechanical model (ϕgeomech). The relationship between fluid porosities and geomechanical porosities would be governed by the following damping equation during the calculations: [00150] Fluid porosity (ψtiuido) can be calculated using Equation 24 and used in fluid flow equations, while geomechanical porosity (ϕgeomec) can be calculated in the geomechanical model and T is a prescribed time constant. After a sudden change in (ϕgeomec), the relative difference between ψfiuido and ϕgeomec may be less than 1%, after a period of 5T. The value of T can be chosen so that the calculations are stable and can be chosen to have a time interval as short as possible, for example, making T have a value that is less than the time period of a time interval in the calculation or the total calculation time. For example, if a calculation is expected to take several days, then T may be in the order of minutes; if a calculation is expected to take a few minutes, then T may be in the order of milliseconds. The T value can be adjusted to about one minute for calculations on many poroplastic materials. In other examples, the value of T can be set to zero (which eliminates all damping). The value of T to be used for a specific calculation will be evident to those skilled in the art. 5.5.2.3 Plastic Flow [00151] The plastic flow, also called flow, denotes a permanent deformation of a material. Once flow into a material begins, the plastic flow may or may not persist. A flow condition is a mathematical representation that marks the transition from elastic to plastic deformation. An assumption in the plastic flow equations is that a single flow condition and a single plastic potential can be used to describe the material of the reservoir and that a single hardening parameter can be used to represent the movement of the flow surface. It can also be assumed that the stress and strain indexes can be expressed as: where the Eijki tensioner indicates the linear elastic properties of the shell material: where À and p are Lame constants. Eq. 26 can be used for an elastic isotropic material. In certain calculations, the plastic model may be limited to isotropic elastic properties. In other calculations, Eqs. 25 and 26 can include other terms that model plastic properties of a material. [00152] It can be assumed that there is a plastic potential G that may be equal to the flow criterion, in some cases, but that may be different from the flow condition for non-associated flow. In addition, it can be assumed that the plastic deformation indices are given by: where the scalar value 2> 0, which is not related to the Lame constant, is related to certain constraints that can be applied to the material properties, that is, the flow condition. The plastic G potential can be used to determine the directions of the plastic deformation indices, while the scalar can be used to determine their magnitudes. [00153] An effective plastic deformation (εp) can be defined as: and the relationship between the rate of change of the effective plastic tension and parameter 2 is: [00154] If a unit tensor ny is defined as: the plastic strain index tensor can be written as [00155] It can be assumed that a simplified form of the flow condition can be written in terms of the stresses (or) and a hardening parameter (K) as: where F is negative in the elastic state and zero in the plastic state and K is a function of plastic deformations (deformation-hardening). During the plastic deformation, the flow condition satisfies the relationship: [00156] Equation 33 can be combined with Equations 25, 26 and 27 to arrive at the expression: [00157] The constitutive equation relates the stress indices to the strain indices when the coefficients depend on the elastic properties, current stresses, current plastic deformations and a hardening parameter. In the case of a type-associative flow calculation, it can be assumed that F = G, which makes the constitutive equation symmetric. If the nonlinear flow equations are solved using an implicitly coupled calculation, the Jacobian for the system of equations can be slightly non-symmetric, even when the constitutive equations are symmetric. 5.5.2.4 Hardening Models [00158] An examination of the requirements for subsequent plastic deformation and the stress-strain ratio of the material can provide an indication of whether or not the plastic flow will persist. A material can ideally be plastic or subject to deformation-hardening. An ideally plastic material (such as, but not limited to, structural steel) exhibits a flow condition that remains unchanged by plastic deformation. However, many materials are altered by inelastic deformation (called deformation-hardening) and the flow condition can be modified as the material's resistance to flow increases. [00159] The geomechanical model includes hardening models that can be used to model the plastic deformation of material within the geomechanical reservoir system. Examples of hardening models include the Drucker-Prager model with shear-hardening or “cap” hardening, the modified Drucker-Prager model with tabular hardening, the modified Bigoni-Piccolroaz model and the Matsuoka-Nakai model. [00160] Drucker-Prager with Shear Hardening [00161] The flow condition for the Drucker-Prager equation with shear hardening and positive tensile stress can be expressed as: where all scalar parameters are non-negative constants, α is a constant, m is an exponent and F is the flow constant. Values for parameters α, k, m and F corresponding to the input parameters can be entered to perform a calculation. The α constant can assume values between 0.0 and 1 / Vã. If the value of α is zero, then the Drucker-Prager model becomes a Von Mises plasticity model. The α parameter is related to the friction angle that is used for a Mohr-Coulomb model. The effective εpplastic deformation is a hardening parameter that can be calculated and communicated, for example, to a user. [00162] After hardening, that is, an increase in εp, the flow surface can be considered to move from an original surface position to a final surface position, as shown in Fig. 11 A. In this calculation, the first stress invariant h (the first total stress invariant) is negative in compression (A = summation in k). [00163] The plastic flow equation for the Drucker-Prager model can be represented by: where εp = λ√3α2 + 1/2 and α in Equation 35 is the same α used in the flow equation, if we assume associated flow. [00164] Drucker-Prager with “cap” Hardening [00165] A Drucker Prager model with “cap” hardening can have two flow surfaces. A drainage surface is the non-compliant Drucker-Prager rupture surface with (perfectly plastic) hardening given by: where all scalar parameters are non-negative constants and tensile stresses are positive, α is a parameter and F is the flow constant. The second flow surface is an elliptical hardening “cap” surface with the formula: [00166] The variable X <RT is a function of the plastic volumetric deformation (compaction is negative) and is the value of the first stress invariant h in which the elliptical “cap” intercepts the axis of hA variable L <0 is also a function of the plastic volumetric deformation and is the value of h at which the elliptical “cap” intersects the Drucker-Prager surface Fs. The elliptical cap is vertical at the intersection with the li axis and horizontal at the intersection with the Drucker-Prager rupture surface. The following equations can be used to relate X to the effective plastic volumetric deformation and the constraint reinforcing the zero slope for the "cap" on the Drucker-Prager runoff surface: [00167] Two “cap” hardening models can be obtained from Equations 39 and 40. The first “cap” hardening model uses the logarithm in Equation 39 for “cap” hardening. The second hardening model “cap” uses the tabular function H () of Equation 39, where H () increases in a way - p strictly monotonic and H (0) = 0. The effective plastic volumetric deformation εve the variables L and X are all zero or negative in Equations 39 and 40. The initial value of εv may be different from zero for a calculation in which the initial magnitude of X exceeds the magnitude of Xo [00168] The movement of the flow surface is shown in Fig. 11B, where the solid line is the original surface and the dotted line is the location of the surface after hardening, and Li and l_2 are the values of L for the surfaces of initial and subsequent flow. [00169] Using the plastic equations, it can be shown that the volumetric plastic strain index can be expressed in terms of parameter 2 from the equation where h <L <0 whenever the plastic deformation occurs in the “cap” and év - ° e 2> 0. Hardening parameters that can be calculated and communicated, for example, to a user, are L, εve Ã. The plastic flow equation on the shear surface can be the same as previously described for the conventional Drucker-Prager model, while the flow on the “Cap” can be given by: [00170] Modified Drucker-Prager with Tabular Hardening [00171] The flow condition for the modified Drucker-Prager equation with shear hardening and positive tensile stress is where α is a constant, K is the deviation, F is the flow constant and 0 is TT / 3 for a uniaxial compression test and can be defined by H () is a tabular function of the effective plastic deformation, h is the first stress invariant and J2 and J3 are the second and third invariants of the total stress diverter: where Sy is the deviation voltage (minus the mean voltage), δy is the Kronecker delta (which is 1 when i = je 0 otherwise), J2 is added to i and je is added to i, je k. The modified Drucker-Prager model is reduced to the Drucker-Prager model if the material parameter K is set to 1.0 and if the table values of H () satisfy Qva | or q0 parameter K can vary from 0.78 to 1.0. The hardening parameter of the effective plastic deformation εpp can be communicated, for example, to a user. The constants α, K and F can be entered. [00172] The effect of parameter K is illustrated in Fig. 12A, which shows graphs of the Drucker-Prager flow surface modified in the octahedral (deviatoric) plane for four K values in the octahedral plane (and for a constant li). In the four graphs in Fig. 12A, it can be assumed that the compression is positive and all have been normalized to cross 1.0 on the positive vertical axis. The conventional Drucker-Prager has K = 1.0 and is circular in the plane. The modified Drucker-Prager surface is no longer convex for K values below 0.78. Certain calculations can be restricted so that K values are greater than or equal to 0.78. [00173] The plastic flow equation for the modified Drucker-Prager equation is where Ty is defined by and Of in the flow equation is the same used in the flow equation when an associated flow is assumed. If it is assumed that Tn = 0 and SÜ = 0, setting Of equal to zero in Equation 49 results in the plastic volumetric deformations being equal to zero for a calculation. [00174] Modified Bigoni-Piccolroaz with Hardening Tabular [00175] The modified Bigoni-Piccolroaz (MBP) model is a modification of the flow criteria in D. Bigoni et al. (2004) “Yield criteria for quasi-brittle frictional materials,” Inti. J. of Structures 41: 2855-2878. The yield equation for the MBP model with positive tensile stress can be expressed as: where α is a constant, I "is the flow constant, y θ θ deviation and where 0 <α <1 / √3, 0 <β <1 Г> 0e0 <ϒ <1 Sg0 material constants and H ( ) is a tabular function of the effective plastic deformation. Parameters J2 and J3 are invariant of the stress diverter and li is the first stress invariant. The constant A can be chosen so that Acos (ψ) is 1.0 for a uniaxial compression test. The MBP model uses deformation-hardening to model plastic hardening and it can be assumed that the flow surface expands evenly (without rotation) as it hardens. The constants α, β, T and y, which can be input parameters, control the shape of the flow surface in the main stress space The angle of the flow cone in the main stress space is controlled by a, while the location of the cone tip on the li axis is given by I7a before The flow surface shape in the octahedral plane (for constant li) can be controlled by parameters β and y. Several common flow surfaces can be reproduced for specific choices of MBP parameters. Common runoff surfaces, which are special cases of MBP runoff are Von Mises, Drucker-Prager, Mohr-Coulomb, Lade, Tresca and Matsuoka-Nakai. [00176] The parameters β and y can be modified to determine the shape of the flow surface in the octahedral plane. The flow surfaces that result in different values of β and y are shown in Fig. 12B. In the graphs of Fig. 12B, it can be assumed that the compression is positive and all graphs have been normalized to cross 1.0 on the positive vertical axis. The MBP flow equation does not exhibit the loss of convexity that can occur for the modified Drucker-Prager flow equation. For example, in the case where β = 0.0 and y = 1.0, the MBP flow equation generates a triangle in the octahedral plane (compared to the non-convex surfaces for the modified Drucker-Prager flow equation) , when K = 0.60). [00177] The plastic flow equation for the MBP model is and at in the flow equation is the same α used in the flow equation assuming an associated flow. If cif is set to zero in Equation 52, the plastic volumetric deformations are zero for the calculations. 5.5.2.5 Sand Production Model [00178] One or more parameters of a hardening model for calculations of plastic deformation can be used for the calculation of the sand production model. In calculations using the sand production model, at least one sand production criterion can be applied to one or more grid cells. It can be assumed that a grid cell fails when a critical value is reached for: (1) the total deformation invariant, or (2) the plastic deformation invariant, or (3) the maximum effective stress at the center of a grid cell. With each of these sand production criteria, the failure of a grid cell represents sand formation from reservoir material (resulting in a cavity). [00179] In the case of criterion (1), the total deformation invariant can be expressed as: [00180] In the case of criterion (2), the plastic deformation invariant can be expressed as: [00181] In certain calculations, the plastic deformation invariant criterion may not be applied if the hardening model is the Drucker-Prager hardening model. [00182] In the case of criterion (3), the maximum effective voltage, the main maximum effective voltage can be calculated and its value can be compared with the input value (where the tensile stresses are assumed to be positive). [00183] Criterion (2), the plastic deformation invariant, can be used to take into account the production of transient sand in a cell before the total failure of that cell. If plastic deformation is critical for total lim cell failure, then the fraction of sand produced from a cell before the total collapse of the cell can be represented by a production function. The production function can be a function f (x), where f (x) = 0, when x = 0, and f (x) = 1, when x = 1, and where x is a function of the critical plastic strain (x = g (εplim)). The function f (x) can be any monotonic function of x, including, but not limited to, a fraction function, a power function, a sine function, a cosine function, a logarithmic function, an exponential function, a sigmoid function or any combination thereof. In some examples, x may be a function of the critical plastic deformation (εplim) and the plastic deformation invariant (εp), such as, but not limited to, an εp / εplim ratio. As an example, the production function can be given by: where m is an exponent. 5.5.2.6 Calculation in the Geomechanical Model [00184] The calculation of the geomechanical model can be in the 3D grid (which can be the grid used for the reservoir model). For example, a conventional finite element method can be used for geomechanical equations where stresses are integrated at the eight Gaussian points inside each element and fluid pressures are integrated at the center of each element. Discretization produces a 27-point model for displacements when eight elements share a common corner and there are three displacements at each node. In certain calculations, displacements are the main variables for geomechanical equations when displacements are evaluated at the corners (based on nodes) of hexahedral elements. 5.5.3 Thermal Model [00185] The thermal model calculates temperature changes that occur when hot or cold fluids are injected into or produced from a reservoir and also calculates the conduction from a drilling well that can circulate hot / cold fluids, but is not , in practice, injecting or producing fluids. Temperature calculations can include porous and geomechanical flow, but can be configured to not include steam injection. It can be assumed that the injected water is in the liquid phase. [00186] The thermal model calculates temperature changes that arise due to conduction and convection in the reservoir, injection / production of hot or cold fluid in wells, conduction in wells and thermal and mechanical interactions between fluids and solids. These mechanisms can be combined into a single energy equation that is solved together with the porous flow equations and solid strain equations. The energy equation can be formulated in terms of a Lagrangian approach to the porous solid and any movement of fluid can be relative to the movement of the porous solid. [00187] With this type of formulation, the mass of the porous solid can be constant when evaluating the energy balance for an element of the reservoir, while the mass of fluid changes when the fluids flow in and out of the porous solid. 5.5.3.1 Combined Accumulation Term [00188] The combined expressions for the change of energy in the fluid and solid (when the geomechanical system is included in the calculation) can be given by: where ψ is the porosity with respect to the apparent undeformed volume, α is the phase, pa is the phase density, Sa is the phase saturation a, ha is the specific enthalpy of phase α, err is the coefficient of volumetric expansion in stress / strain equations, K is the permeability, given the volumetric expansion coefficient in the porosity equation and Toé the fluid temperature. Conduction and convection terms can affect the combined accumulation term. In Equation 56, an approximation can be made in which the phase pressures are equal, the term P can be approximated by Pe the term Cr can be -Cv + avMTo & heat capacity Cv is the heat capacity for a porous solid measured in an interstitial volume constant and constant apparent volume. It can be shown that the heat capacity Cr, as defined above, is the heat capacity at a constant fluid pressure and constant apparent volume. Equation 56 is the general formula for the total index of internal energy change in an element and is used for the accumulation term when terms including the temperature of fluid are considered in the calculation and the calculation includes the thermal and geomechanical models. [00189] Equation 56 takes on simpler forms if the fluid temperature To is considered in the calculation. [00190] In the case of calculations that include thermal and geomechanical models and terms including fluid temperature To and Vo / 7 are not considered in the calculation, the simplified expression is [00191] In the case of calculations that include the thermal model, but do not include the geomechanical model, the accumulation term is: while for calculations that include the thermal model, but do not include the geomechanical model and the internal fluid energy is approximated by the fluid enthalpy, the expression is: 5.5.3.2 Energy Equation [00192] The choice of the accumulation term above determines the equations that are included in the final energy conservation equation. The energy equation can be expressed as where a constant temperature boundary condition can be applied to injection wells. Since Equation 60 does not contain any mechanical terms, with the exception of transport, temperatures in the grid blocks should not be affected by the expansion or compression of a fluid or solid phase when enthalpies are just temperature functions. However, Equation 60 still takes into account the heat of vaporization from a liquid to vapor because latent heat can be included in the difference between the enthalpy of a component in a liquid phase and its corresponding enthalpy in the vapor phase. [00193] If a simulation includes geomechanical calculations and terms including the fluid temperature To are considered in the calculation, then the general energy equation can be given by: [00194] The term vestá present on both sides of the energy equation, that is, in the accumulation term and also as work done on the set of solid particles and, consequently, σtiεv is not included in Equation 61. 5.5.3.3 Modifications of Fluid Densities and Viscosities [00195] There are several options for calculating densities and viscosities for temperature-dependent fluids and the methods available for calculating densities and viscosities vary depending on the phase behavior model and the numerical technique being used. In certain examples, fluid properties can be calculated directly as a function of temperature and pressure. In certain other examples, fluid properties can be calculated as a function of current cell pressures and initial cell temperatures (isothermal flash) and modification functions can be used to correct differences between initial cell temperatures and current cell temperatures. cells. [00196] Fluid properties can be calculated directly as a function of temperature and pressure for implicit single-phase executions, implicit two-phase executions and implicit compositional executions. These calculations involving the geomechanical model can be coupled iteratively or fully coupled. [00197] The viscosities and densities that are calculated in black oil, K value and compositional PVT packages can be modified to take into account changes in fluid properties due to temperature changes. The initial reservoir temperature and the current fluid pressure during a calculation can be used to calculate the fluid properties in each of the PVT (isothermal flash) packages, and then these properties are modified by multiplying them by modification factors provided by the user, ie P (P, TO) P (T em which can be inserted as a table and XPX) can be calculated in the PVT package. [00198] In the case of many thermal-geomechanical studies, the interest can be focused mainly on how the temperature affects the fluid viscosity and thermal stresses. For these types of applications, modification factors or a direct calculation of fluid properties can be used and the results can be very similar in both techniques. However, in the case of applications where temperature has a strong effect on density and / or phase separation, then the direct calculation of fluid properties can be used as a function of temperature and pressure. 5.6 Interrelationship of the Various Models [00199] A system and method can couple the interactions between the different models of geomechanical reservoir systems, for example, between porous flow, displacements and geomechanical stresses, sand production and cavity formation. For example, the geomechanical solution can influence the porous flow calculations through the terms of porosity and permeability. In another example, the fluid solution in the reservoir can affect geomechanical calculations through its role in effective stresses. In yet another example, the fluid solution can affect the geomechanical calculations involving the sand production model through normal traction that occurs on the face of the cavity and can affect the flow in the reservoir. 5.7 Implementation of Sand Production Prediction 5.7.1 Calculations Including the Sand Production Model [00200] Fig. 13 illustrates several considerations in an exemplary calculation of the prediction of sand production in a time interval of a calculation in which a cavity was generated between the reservoir cover soil (rock covering the area of interest in the well) and the underlying soil (rock underlying the area of interest in the well). At the end of each time interval, a check can be made for failed grid cells and a cell can be removed immediately from strain calculations (for example, a calculation involving a geomechanical model) as soon as it is identified like a failed cell. A failed cell is a cell for which a sand production criterion has been met (see Section 5.5.2.5 above), which can be used to indicate the point where there is a material failure, that is, the material becomes granulated without binder to form sand and generate a cavity. After a cell fails and becomes part of the cavity, it can be excluded from the stiffness matrix to calculate the displacements. As shown in Fig. 13, it can be checked whether the flawless reservoir cells in front of the interface between an existing cavity and flawless cells have flaws at each time interval. As also illustrated in Fig. 13, it can be assumed that the tensile condition at the interface between an existing cavity and flawless cells can be equal to the fluid pressure in the cavity cells at the interface. The porous flow models and conventional geomechanical models can be applied to intact grid cells without fail. [00201] The volume of sand in the failed cell can then be added to the volume of sand production reported to the nearest production well. In some calculations, the volume of sand in a cell can be the initial apparent volume of that cell. The failed cell can be retained in the fluid flow calculations to allow flow between the generated cavity wall and the well, and the cell porosity can be kept constant. In some calculations, permeability can be increased in the cavity to minimize the pressure drop between the well and the cavity wall (for fluid flow). After the cell is removed from the deformation calculations, its properties no longer affect the deformations, except that the cells with failure immediately adjacent to the cavity can exert a compression load on the porous solid without fail (that is, a cell without failure) and this compressive load on the cavity wall can be directly related to the fluid pressures in the failed cells adjacent to the cavity wall. [00202] Cavity generation calculations can be physically and / or numerically unstable. To avoid physical instability, some hardening can be added to the end of a curve that is almost elastic-perfectly plastic and a time interval can be defined for the growth of an unstable cavity (ie the growth rate for the cavity can be be directly related to the time scale set, such as, by way of non-limiting example, one minute). In order to minimize numerical instabilities, several parameters of the models can be adjusted to certain values. For example, the time constant T of Equation 24 (which relates changes in fluid porosities with changes in geomechanical porosities) can be adjusted to a value other than zero to minimize instabilities, since a value of zero removes all damping. A higher value of the time constant T can improve numerical stability. In another example, the permeability attributed to grid cells that fail (indicating sand formation) can be adjusted to a value that minimizes numerical instability. As a non-limiting example, permeability can be adjusted to a value on the order of about 1000 Darcies. Lower or higher values can be set for permeability to increase numerical stability. 5.7.2 Solution of the System of Partial Differential Equations [00203] The Jacobian for the complete system of equations for the models included in the calculation can be completely expanded and all equations can be solved simultaneously in the linear resolution agent. The various mechanisms can be coupled using techniques unidirectionally coupled, explicitly, iteratively or totally. The most stable coupling technique can be used for predictions of sand production. [00204] The equations discussed in the previous sections can be implemented in a single program and the program calculates an implicit solution for all equations, that is, a regressive Euler technique can be used to approximate the time derivatives and all variables and coefficients equations can be evaluated at the end of the time interval. Non-limiting examples of primary variables are fluid pressures in the reservoir, displacements in the reservoir and the boundary of the reservoir (if not restricted), fluid pressures in a fracture (or fluid volumes in the fracture if partially saturated), the pressure of well and at least one parameter of the production function. The system of equations can be solved using a Newton-Raphson technique where the fully expanded Jacobian is generated for the entire system of equations and incremental corrections are found using a linear resolution agent that includes all solution variables. [00205] The cavity generation can be simulated by calculating flow to a failed cell and calculating the effective stresses and displacements that arise in the adjacent reservoir. Criteria for sand production are then combined with the stress and displacement states and traction condition at the cavity interface to decide whether the cavity has progressed. The cavity can progress through multiple grid cells over a single period of time. Each time slot can start with the cavity setting, geomechanical state, fluid state (if applicable) and thermal state (if applicable) of the previous time slot. The calculation then iterates for a converged solution for the conservation of mass and tension balance equations assuming that the cavity does not propagate. After obtaining a new convergent solution, the sand production criteria can be checked to see if the critical value of the sand production criterion has been reached in any cells. These cells can be treated as described in Section 5.7.1. This sequence of iterations can be continued in the equations and checking sand production criteria for a series of times or until a continuation of the cavity progression over a period of time is not seen before moving on to the next period of time. 5.8 Examples of Appliances and Computer Program Implementations [00206] The methods disclosed herein can be implemented using an apparatus, for example, a computer system, such as the computer system described in this section, according to the following programs and methods. This computer system can also store and manipulate data indicating physical properties associated with a geomechanical reservoir system, the Jacobian fully expanded to the complete system of equations for the models included in the calculation, the Jacobian solution fully expanded, the predictions generated fractures or measurements that can be used by a computer system implemented with the analytical methods described here. The systems and methods can be implemented in various types of computer architectures, such as, for example, in an individual general purpose computer or in a parallel processing computer system, or in a workstation, or in a networked system (for example, a client-server configuration, as shown in Fig. 14). [00207] As shown in Fig. 14, the modeling computer system to implement one or more methods and systems disclosed herein can be connected to a network connection that can be, for example, part of a local area network (“ LAN ”) to other local computer systems and / or part of a wide area network (“ WAN ”), such as the Internet, which is connected to other remote computer systems. [00208] The modeling system comprises any methods as described here. For example, a software component can include programs that cause one or more processors to implement steps that consist of accepting a plurality of parameters indicating physical properties associated with the geomechanical reservoir system and / or the generated fracture predictions, and storing the parameters indicative of the physical properties associated with the geomechanical reservoir system and / or the fracture predictions generated in memory. For example, the system can accept commands to receive parameters indicative of physical properties associated with the geomechanical reservoir system and / or parameters of a generated fracture prediction, which are entered manually by a user (for example, through the user interface. ). The programs can cause the system to extract parameters indicative of physical properties associated with the geomechanical reservoir system and / or parameters of a fracture prediction generated from a data storage medium (for example, a database). This data storage medium can be stored in a mass memory (for example, a hard disk drive) or other computer-readable medium and inserted into the computer's memory, or the data storage medium can be accessed by the computer system through the network. 6. RESULTS 6.1 Sand Production Prediction Tests [00209] Sand formation tests were carried out on samples of saturated cores of salt water and kerosene. The prediction of sand production was applied to the results of sand formation tests. The effects of temperature were not considered. It was assumed that single-phase simulations could be used to simulate the reservoir system. 6.2 Data Collection [00210] The geometry of the grid, information on the material behavior, boundary conditions and permeability function are discussed below. 6.2.1 Geometry / Grid [00211] The tested samples were approximately 100 mm long and had a circular diameter of approximately 100 mm. The central hole had a diameter of approximately 20 mm. The exact dimensions of the salt water and kerosene test samples are given in Table 1. Fig. 15 shows the core samples used for the salt water test, including the sample dimensions. During the test, press plates were used to apply the axial load. A hollow cylinder test was applied (as shown in Fig. 3). [00212] An axisymmetric grid of cells was configured for the calculations, based on the dimensions indicated in Table 1. The rubber sleeve was not represented in the cells of the calculation grid. However, the simulation model was configured to include the two plates, whose radial dimensions correspond to those of the sample and whose axial dimension was 20 mm. [00213] Table 1: Sample Dimensions. 6.2.2 Material behavior [00214] Physical constants of the materials were determined using data from a triaxial assay. The material response measured and calculated using the models are in close agreement both in the case of cores saturated with salt water (shown in Fig. 16) and in the case of cores saturated in kerosene (shown in Fig. 17). The physical constants used to model the sample saturated in salt water are given in Table 2. [00215] Table 2. Physical properties for the sample saturated in salt water. 6.2.3 Boundary conditions [00216] The simulations were performed by solving the fully coupled fluid flow and geomechanical models and, therefore, the boundary conditions for both fluid flow and geomechanics had to be defined (shown in Fig. 18). In the case of geomechanical boundary conditions, a constant stress was attributed to the top, outer and inner surface (boundary) of the hollow cylinder (indicated by arrows on the panel (A) of Fig. 18). The geomechanical model was limited in terms of displacement at the bottom to prevent the numerical model from moving and rotating (that is, the entire surface is constrained to not move in the z direction). To simulate the fluid flow test conditions (see panel (B) of Fig. 18), a pressure drop was applied to the hollow cylinder, the pressure on the outside of the cylinder being greater than the pressure inside the cylinder. . No fluid flow boundary conditions were applied to the top and bottom of the cylinder. The porosity of the spacers was set to zero, which, in effect, also makes the permeability of the spacers equal to zero. [00217] The pressure inside the hollow cylinder, as well as the stresses on the inner surface of the hollow cylinder, were kept equal to 0 psi. The simulation used absolute pressure values, therefore, it was defined that a zero pressure value was 14.7 psi to take atmospheric pressure into account. Both the stresses outside and the pressures outside the cylinder varied depending on the flow and stress regime applied in the test performed on the sample saturated in salt water. Fig. 19 shows the confinement stress, flow and flow pressure measured and subjected to discretization vs. time for the saturated salt water sample. 6.2.4 Permeability function [00218] Using the test measurements, the permeability was calculated for a homogeneous radial isotropic geometry. Darcy's law for radial flow to the core sample geometry can be expressed as: where k is the permeability, Q is the flow rate, p is the viscosity, right is the inner radius, r is the outer radius, h is the height and ΔP is the pressure drop. The permeability was determined using the dimensions of the core samples and Equation 62. The permeability distribution resulting from this calculation showed that the permeability was not constant as a function of flow (as shown in Fig. 20). Specifically, as the flow increases, the permeability decreases (it was observed that there was an influence on the quality of pressure measurements). An average permeability was determined by adapting a power law curve to a graph of the permeability calculated as a function of the confinement stress (Fig. 21). This power law function was used to calculate the pressure to achieve the flow. In other words, the pressures were changed so that, given the permeability function derived from the adaptation of the power law, the desired flow rate was obtained. [00219] Generally, permeability is a function of cell properties and not boundary conditions, such as confinement stress. The assumption made was that the mean stress can be used as a representation of the confinement stress, the mean stress being calculated using: where the mean is the average stress and σi to 03 are the main stress. This assumption meant that the calculated flow and the measured flow were not identical. The flow measured in the tests was compared with the flow predicted by the numerical simulation. A correction factor was applied to the external boundary pressures to make the flows coincide. The calculated pressures were multiplied by 1.035 to obtain the coincidence of flows (shown in FIG. 22), which meant that only a correction factor of 3.5% was used. This coincidence of flows was obtained for a sample without sand production. Flow rates may vary depending on the volume of sand produced, as discussed in Section 6.3 below. 6.3 Analysis [00220] The calculated sand coincidence and its sensitivity to certain parameters are discussed below. A question is also addressed as to whether a sustainable and stable cavity can be created. 6.3.1 Sand Production Model Calculation [00221] The following inputs for the sand production model were calculated: 1. value of a critical plastic deformation - the critical value of an effective plastic deformation that causes the rock to fail. 2. exponent - exponent of the production function used to calculate the amount of sand production from a cell before reaching the critical plastic deformation value. 3. hardening model parameters - the triaxial test data were matched (adapted), using a modified Bigoni-Piccolroaz (PAM) version of the hardening model. [00222] To adapt the MBP hardening model to the triaxial test data, in addition to using the parameters given in Table 2, it is possible to make the MBP model allow a user to define the shape of the flow surface in the octahedral plane (constant h) (see Figs. 23A and 23B). The shape of the flow surface in the octahedral plane is defined by two parameters for the MBP hardening model: deviation (y) and beta (β). Both the deviation and beta are constant having a value between 0 and 1.0. Setting 0 for beta and 0 for deviation in the MBP model produces a flow surface in the octahedral plane of the material model Drucker Prager (shown in Fig. 23A). The definition of 0 for beta and 0.95 for the deviation in the MBP model produces a flow surface in the octahedral plane of the Type Lade material model (shown in Fig. 23B). [00223] Flow occurs when a sample reaches the flow surface. Thus, a linear elastic behavior is dominant in the inner areas of the shapes shown in Figs. 22A and 22B. From Figs. 22A and 22B, it was found that the probability of a voltage path reaching the flow surface in the deviatoric plane for the Lade Type model is greater than for the Drucker-Prager model. It was expected that more flow would occur with higher deviation values when applying the two material models to the general stress states. It was noted that the flow surfaces in Figs. 23A and 23B produced the same results when applied to triaxial compression tests. 6.3.1.1 Definition of Grid Cells [00224] The calculations were performed on 2D axisymmetric grids. The grid cell size to represent the sample was varied to investigate grid size dependencies on the model results. The sample cell size was 1.284 mm in the vertical direction and varied from 0.05 mm to 0.01 mm to 0.005 mm in the horizontal direction. The measurements were accurate with a tolerance of 0.1 cm3, so small cell sizes were used in the model to capture this level of accuracy. The cell size in the spacers was 2 mm in the vertical direction. This resulted in models with 80300, 401100 and 802100 cells, respectively. Run times progressively increased from about 2 hours to ~ 14 hours to approximately 30 hours. Fig. 24 shows the progression of the model's results as a function of cell size. No visible changes were observed in the predicted sand production volume from the grid with a horizontal cell size of 0.01 mm to 0.05 mm. Run times were reduced below 14 hours, an additional grid cell model was created that had the same vertical cell dimensions as the previous ones, however, only the first 7 mm of the sample were cross-linked with a size of 0.01 mm cell. In addition to the initial 7 mm, the cell size was increased to 0.11035 mm. This graduated grid cell model had 100,000 cells and ran in approximately 3 ~ 4 hours. The calculation with this grid cell model was compared with the results of other grid cell models and it was found that the results for this graduated grid cell model overlapped with those of the thinner grid cell models, indicating that the results are independent of the cell size (Fig. 24). This graduated grid configuration was used for all subsequent runs. 6.3.1.2 Adaptation to the Production Function [00225] The production function used for the adaptation was f = (ε • / inw) m and the simulations were performed by solving the system of partial differential equations of the coupled geomechanical model (including the production function) and the fluid flow model. The best adaptation was obtained for a deviation of 0.5, a beta value of 0, an exponent value (m) of 2 and a critical strain limit of 0.017 (see Fig. 25). The data were adapted for both the massive sand formation part and the gradual sand formation part of the curve. The beginning of the plastic deformation is indicated in Fig. 25. The water flow correspondence is indicated in Fig. 26. Approximate correspondences between the sand formation volume and the water flow were observed. A slight mismatch in the water flow is observed at later times, due to the high permeability values attributed to the cells that were filled with sand. The test results indicated that the increased permeability slightly overestimated the water flow in the numerical model. In addition, it was observed that the critical deformation value of 0.017 indicated that there was sand formation far beyond the initial flow of the sample (Fig. 27), indicating that when the initial flow point is considered as the beginning of formation of sand, the trends of sand formation in a sample can be significantly overestimated. [00226] To obtain the adaptation to the volume of sand produced, the deviation was varied first, with the exponent (m) kept constant with a value of 1. For each deviation value, the critical deformation value was chosen to resemble the runoff value during massive sand formation. To obtain this value, the adaptation was performed using a geomechanical model that contained the hardening model, but not the sand production model. The multiple values of deviation, critical deformation and exponent were used to investigate the unusualness of the solution. [00227] The deviation value was, in the first place, varied from 0.92 to 0.5. As indicated before, the value of critical plastic deformation was determined for each deviation by performing the adaptation using a geomechanical model that contained the hardening model, but not the sand production model. The flow value at the beginning of the massive sand formation was used as the initial value of the critical plastic deformation. Adjustments to the value of the critical plastic deformation were sometimes made to obtain the best possible match given the value of the deviation. A maximum deviation value of 0.92 was used. This value was determined from the Matsuoka-Nakai material model (a version of the MBP model). As shown in Fig. 28, a deviation of 0.5 gave the closest match and, therefore, this value was used in all subsequent simulations. [00228] The variation of the exponent (m) showed that the value of the exponent only influenced the amount of sand produced, but did not change the point at which massive sand formation occurred. The higher values of the exponent reduced the total volume of sand produced (Fig. 29). The variation in the value of the critical plastic deformation showed that this constant governed the beginning of the massive formation of sand, in which higher values of the exponent resulted in the formation of sand in later moments (that is, in higher confinement stresses) (Fig. 30) . The reduction in the value of critical plastic deformation resulted in the occurrence of sand formation at lower confinement stresses. [00229] Fig. 31 shows the evolution of the shape of the material failure from sand formation after the hollow cylinder test has been carried out. Each panel in Fig. 31 is a horizontal strip taken at different depths in the core sample. From Fig. 31, it was observed that the formation of sand was more massive in the core core of the sample (that is, the cavity is deeper in the core of the core) and less at the ends of the core sample. This form was also observed in the numerical simulation runs. FIG. 32 shows the simulated evolution of the material failure shape over time calculated for the numerical simulation (each track in Fig. 31 is captured vertically). Although, initially, the fault was located at the spacer-rock interface, the cavity progressed and the deepest cavity occurred in the core (see Fig. 32), that is, the cavity became deeper as it progressed towards the center of the simulated core, which corresponds to what was observed in the results of hollow cylinder tests. 6.3.2 Cavity Growth [00230] As can be seen from Fig. 32, at 116675 seconds, there was still sand production. To determine whether sand production stabilizes, loading at 116675 seconds (ie, a confinement stress of 46 MPa and a pressure differential of 0.1643398 MPa) was kept constant for an additional 116675 seconds. However, the simulation was stopped if more than 50% of the sample was produced in the form of sand. Fig. 33A shows the sand production (cm3) vs. time (s) and Fig. 33B shows the cavity size for the numerical simulation, with the confinement stress constant at 46 MPa and the pressure differential constant at 0.1643398 MPa from 116675 seconds forward. A continuous grid with a size of 0.01 mm was used. In Fig. B, the gray area indicates the location of the cavity and the dark area indicates intact rock (or spacer). At 126618.836995 seconds, 50% of the sample was produced in the form of sand and the simulation was terminated (see Figs. 33A and 33B). In the volume calculation, the volume of the spacers was also counted. As shown in Fig. 33B, the cavity was very deep in the center and almost completely penetrated the complete simulated rock sample. Thereafter, it was concluded that a confinement stress of 46 MPa and a pressure differential of 0.1643398 MPa resulted in an unstable sand formation. [00231] To determine whether a stable cavity could be developed and sustained, the loading was changed. The material constants used were identical to those in Fig. 25, that is, deviation = 0.5, beta = 0, critical plastic strain value = 0.017 and the exponent (m) = 2. A loading process similar to the one used described in association with Fig. 32. Specifically, the loading reached at 113230 seconds, that is, a confinement tension of 44 MPa and a pressure differential of 1.2555317 MPa, was kept constant for a further 120120 seconds (simulation time max 233350 s). This simulation showed that a maximum of 10.2966 cm3 of sand was produced and that the cavity stabilized (Fig. 34). This observation indicated that there were conditions in which sand production was foreseeable to some extent, but which should not lead to catastrophic failure (due to massive sand formation). 6.4 Summary [00232] The results showed that both the beginning of the sand formation and the massive sand formation of the hollow cylinder sand formation test could be uniquely adapted using a geomechanical model including a production function. The parameters used to adapt the data are deviation, beta, critical deformation and the exponent. For the Landana saltwater test, it was found that the following parameter values resulted in a better match: • Deviation = 0.5 • Beta = 0 • Critical strain value = 0.017 • Exponent = 2 [00233] The results also indicated that there were conditions in which some amount of sand production was expected, but which should not lead to a catastrophic failure (due to massive sand formation). The workflow for the adaptation process in Section 6.3.1 was as follows: 1) Based on the total volume of sand and the precision with which the volume of sand was measured, an appropriate grid size was chosen. For this study, if a total amount of approximately 6 cc of sand was produced, a cell size of 0.01 mm was used to obtain the desired precision. 2) Multiple models with different deviation values were created using the material constants obtained from adaptations to the data of the triaxial test. In the calculations to adapt the test data, the geomechanical model did not include the sand production model. 3) The flow value was recorded in the massive formation of solid sand using the results of step 2. This flow value was defined in order to have a value equal to the critical plastic deformation limit and the exponent value was kept equal to 1. In this step, several passes could be performed with a range of values of critical plastic deformation around the obtained value. 4) Using the results of step 3, the simulation that adapts to the beginning of the massive sand formation was used to vary the exponent value, where an increase in the exponent value was observed to reduce the amount of sand produced. The simulation that best suits the test data was chosen as representative. 7. REFERENCES CITED [00234] All references cited herein are incorporated by reference, in their entirety and for all purposes, in the same way as if each individual publication or patent or patent application was specifically and individually indicated to be incorporated / o here for reference in its entirety and for all purposes. The discussion or quotation of a reference here will not be construed as an admission that such reference is a prior art to the present invention. 8. MODIFICATIONS [00235] Many modifications and variations of this invention can be made without deviating from its spirit and scope, as will be evident to those skilled in the art. The specific examples described herein are provided by way of example only and the invention is limited only by the terms of the appended claims in conjunction with the full scope of equivalents to which such claims are entitled. [00236] As an illustration of the wide scope of application of the systems and methods described herein, the systems and methods described herein can be implemented in many different types of processing devices by program code comprising program instructions that are executable by the processing subsystem. device. The software program instructions can include source code, object code, machine code or any other stored data that can be handled to cause a processing system to perform the methods and operations described here. However, other implementations can also be used, such as firmware or even hardware designed properly configured to execute the methods and systems described here. [00237] The data of the systems and methods (for example, associations, mappings, data entry, data output, intermediate data results, final data results, etc.) can be stored and implemented in one or more different types of computer-implemented data stores, such as different types of storage devices and programming constructs (eg, RAM, ROM, Flash memory, flat files, databases, programming data structures, programming variables, instruction constructs IF-THEN (or similar type), etc.). It should be noted that data structures describe formats for use in organizing and storing data in databases, programs, memory or other computer-readable media for use by a computer program. [00238] Systems and methods can be provided in many different types of computer-readable media including computer storage mechanisms (eg, CD-ROM, floppy disk, RAM, flash memory, computer hard disk drive, etc.) which contain instructions (for example, software) for use in the execution by a processor to perform the operations of the methods and implement the systems described herein. [00239] The computer components, software modules, functions, data stores and data structures described here can be connected directly or indirectly to each other in order to allow the data flow necessary for its operations. It should also be noted that a module or processor includes, but is not limited to, a code unit that performs a software operation and can be implemented, for example, as a subroutine code unit or as a code unit software function, either as an object (as in an object-oriented paradigm) or as a mini-application, or in a computer scripting language or as another type of computer code. Software components and / or functionality can be located on a single computer or distributed across multiple computers, depending on the situation.
权利要求:
Claims (12) [0001] 1. Method for use in forecasting sand production from a geomechanical reservoir system, characterized by the fact that it comprises: receiving (100) data indicative of plastic deformation, gradual sand formation and massive sand formation of the material within the reservoir system geomechanical; compute (102), in a computer system, a value of one or more hardening parameters based on a first adjustment of a hardening model for the received plastic deformation data; wherein said hardening model models the behavior of the plastic deformation of said material within the geomechanical reservoir system; compute (104), in a computer system, a value of a critical plastic deformation based on a second adjustment of the geomechanical model comprising said hardening model for the data received from massive sand formation; compute (106), in a computer system, a value of at least one parameter of a production function based on a third adjustment of the geomechanical model comprising said production function for the received sand formation data received and using a value at least one value of said hardening parameters and said value of said critical plastic deformation; wherein said production function provides for sand production from said geomechanical reservoir system; and sending (108) to a user interface device, a monitor, a computer-readable storage medium, a local computer or a computer that is part of a network, at least one value within said one or more hardening parameters , said value of said critical plastic deformation or said value of at least one parameter of said production function. [0002] 2. Method, according to claim 1, characterized by the fact that said production function is a function f (x), where f (x) = 0 when x = 0 and f (x) = 1 when x = 1 ; and where x is a function of the critical plastic deformation. [0003] 3. Method, according to claim 1, characterized by the fact that said value of at least one parameter of said production function is a value of an exponent of said production function. [0004] 4. Method, according to claim 3, characterized by the fact that the said function of sand production is given by the expression: [0005] 5. Method, according to claim 1, characterized by the fact that said hardening model is a modified Bigoni-Piccolroaz model. [0006] 6. Method, according to claim 1, characterized by the fact that the data of said plastic deformation are obtained from a triaxial test; or the data of said gradual sand formation and the mass sand formation data are obtained from a hollow cylinder test. [0007] 7. Method, according to claim 1, characterized by the fact that said first adjustment of the hardening model for the received data of plastic deformation is obtained by a regression. [0008] 8. Method, according to claim 1, characterized by the fact that: the said second adjustment of said geomechanical model comprising said hardening model for the received data of massive sand formation is obtained by solving a system of partial differential equations that model the geomechanical reservoir system; the system of partial differential equations comprises a reservoir flow model and said geomechanical model comprising said hardening model; in which the system of partial differential equations is coupled through a fully expanded Jacobian equation; and the resolution of the system of partial differential equations includes resolving simultaneously in a single time step the Jacobian fully expanded based on the data received from massive sand formation. [0009] 9. Method, according to claim 1, characterized by the fact that: said third adjustment of said geomechanical model comprising said production function for the received gradual sand formation data is obtained by solving a system of partial differential equations that model the geomechanical reservoir system; the system of partial differential equations comprises a reservoir flow model and said geomechanical model comprising said hardening model; in which the system of partial differential equations is coupled through a fully expanded Jacobian equation; and the resolution of the system of partial differential equations includes resolving simultaneously in a single time step the Jacobian fully expanded based on the gradual sand formation data received. [0010] 10. Method, according to claim 1, characterized by the fact that at least one value sent among said one or more hardening parameters, said value of said critical plastic deformation or said value of at least one parameter of said production function, is visually displayed on the user interface device or on the monitor. [0011] 11. Method for producing sand from a geomechanical reservoir system, characterized by the fact that it comprises operating said geomechanical reservoir system according to a result of the method as defined in claim 1. [0012] 12. Method for operating a geomechanical reservoir system to control sand production from the geomechanical reservoir system, characterized by the fact that it comprises: computing a value of at least one parameter of a production function according to the method as defined in any one of claims 1 to 10, the at least one parameter value of a production function being an operating parameter and the value of at least one operating parameter indicating a condition for the production of stabilized sand from the geomechanical reservoir system; and operating said geomechanical reservoir system according to the value of at least one operating parameter.
类似技术:
公开号 | 公开日 | 专利标题 BR112012006055B1|2020-09-15|METHODS FOR USE IN THE PROVISION OF SAND PRODUCTION OF A GEOMECHANICAL RESERVOIR SYSTEM, FOR SAND PRODUCTION OF A GEOMECHANICAL RESERVOIR SYSTEM, AND TO OPERATE A GEOMECHANICAL RESERVOIR SYSTEM TO CONTROL THE SAND PRODUCTION SYSTEM OF THE GEOMECHANICAL SYSTEM OF THE SAND SYSTEM OF THE SYSTEMATIC SAND SYSTEM. US9747393B2|2017-08-29|Methods and systems for upscaling mechanical properties of geomaterials US9164194B2|2015-10-20|Method for modeling deformation in subsurface strata Kim et al.2013|Development of the T+ M coupled flow–geomechanical simulator to describe fracture propagation and coupled flow–thermal–geomechanical processes in tight/shale gas systems Hui et al.2013|The upscaling of discrete fracture models for faster, coarse-scale simulations of IOR and EOR processes for fractured reservoirs WO2010047859A1|2010-04-29|Method for modeling deformation in subsurface strata Lin et al.2017|A criterion for evaluating the efficiency of water injection in oil sand reservoirs Alonso et al.2013|Joints in unsaturated rocks: Thermo-hydro-mechanical formulation and constitutive behaviour Bai et al.2016|Hydraulic Fracture Modeling Workflow and Toolkits for Well Completion Optimization in Unconventionals Xiong et al.2014|Modification of thermo-elasto-viscoplastic model for soft rock and its application to THM analysis of heating tests Li et al.2021|Stress-dependent fracture permeability measurements and implications for shale gas production Deng et al.2019|Influence of sand production in an unconsolidated sandstone reservoir in a deepwater gas field Zhang et al.2019|Inference of in situ stress from thermoporoelastic borehole breakouts based on artificial neural network Wu et al.2021|Advances in statistical mechanics of rock masses and its engineering applications Huang et al.2017|Poro-viscoelastic modeling of production from shale gas reservoir: An adaptive dual permeability model Jia et al.2020|A nonlinear elasto‐viscoplastic model for clayed rock and its application to stability analysis of nuclear waste repository CN108572401B|2020-04-03|Construction method of fracture-cavity combined model and method for detecting deformation of reservoir fracture-cavity Dubinya et al.2015|Two-way coupled geomechanical analysis of naturally fractured oil reservoir's behavior using finite element method Shin et al.2015|On computation of strain‐dependent permeability of rocks and rock‐like porous media Wang et al.2020|Study of gas emission law at the heading face in a coal‐mine tunnel based on the Lattice Boltzmann method Zheng et al.2014|Investigation of Coupled processes and impact of high temperature limits in argillite rock Li et al.2021|Modeling hydraulic fracture propagation in a saturated porous rock media based on EPHF method Tanaka et al.2013|A novel approach for incorporation of capillarity and gravity into streamline simulation using orthogonal projection Li et al.2017|Investigation of wellbore breakouts in deviated wells-A 3D numerical modeling approach Wu et al.2020|Modeling Thermal-Hydrologic-Mechanical Processes for EGS Collab Thermal Circulation Tests using Embedded Discrete Fracture Model
同族专利:
公开号 | 公开日 EP2478457B1|2019-06-05| WO2011035146A2|2011-03-24| AU2010295468A1|2012-04-12| US20140122035A1|2014-05-01| US9026419B2|2015-05-05| EP2478457A4|2017-03-22| EA201270432A1|2013-01-30| WO2011035146A3|2011-08-25| CN102575510A|2012-07-11| US20110061860A1|2011-03-17| AU2010295468B2|2015-07-02| EP2478457A2|2012-07-25| CA2774045A1|2011-03-24| BR112012006055A2|2016-03-29| CN102575510B|2015-01-28| US8548783B2|2013-10-01|
引用文献:
公开号 | 申请日 | 公开日 | 申请人 | 专利标题 US6980940B1|2000-02-22|2005-12-27|Schlumberger Technology Corp.|Intergrated reservoir optimization| FR2869116B1|2004-04-14|2006-06-09|Inst Francais Du Petrole|METHOD FOR CONSTRUCTING A GEOMECHANICAL MODEL OF A SUBTERRANEAN ZONE FOR TORQUE TO A RESERVOIR MODEL| US20080167849A1|2004-06-07|2008-07-10|Brigham Young University|Reservoir Simulation| AU2005259253B2|2004-06-25|2008-09-18|Shell Internationale Research Maatschappij B.V.|Closed loop control system for controlling production of hydrocarbon fluid from an underground formation| US8301425B2|2005-07-27|2012-10-30|Exxonmobil Upstream Research Company|Well modeling associated with extraction of hydrocarbons from subsurface formations| US7660670B2|2006-10-27|2010-02-09|Schlumberger Technology Corporation|Sanding advisor| US7653488B2|2007-08-23|2010-01-26|Schlumberger Technology Corporation|Determination of point of sand production initiation in wellbores using residual deformation characteristics and real time monitoring of sand production| US20090125280A1|2007-11-13|2009-05-14|Halliburton Energy Services, Inc.|Methods for geomechanical fracture modeling|US8374836B2|2008-11-12|2013-02-12|Geoscape Analytics, Inc.|Methods and systems for constructing and using a subterranean geomechanics model spanning local to zonal scale in complex geological environments| WO2012051184A2|2010-10-14|2012-04-19|Baker Hughes Incorporated|Predicting downhole formation volumetric sand production using grain-scale rock models| US20120203524A1|2011-02-09|2012-08-09|Conocophillips Company|Quantitative method of determining safe steam injection pressure for enhanced oil recovery operations| WO2012118867A2|2011-02-28|2012-09-07|Schlumberger Technology Corporation|Method to determine representative element areas and volumes in porous media| CA2858763C|2011-12-20|2020-06-16|Shell Internationale Research Maatschappij B.V.|Automated calibration of a stratigraphic forward modellingtool using a neighborhood algorithm with explicit escape clauses| US10077639B2|2012-06-15|2018-09-18|Landmark Graphics Corporation|Methods and systems for non-physical attribute management in reservoir simulation| AU2013315927B2|2012-09-13|2017-09-21|Chevron U.S.A. Inc.|System and method for performing simultaneous petrophysical analysis of composition and texture of rock formations| EP2923225A4|2012-11-20|2016-10-12|Stochastic Simulation Ltd|Method and system for characterising subsurface reservoirs| CN103077556B|2013-02-04|2016-07-06|重庆大学|The Three-dimension Numerical Model method for designing of sand production| US9189576B2|2013-03-13|2015-11-17|Halliburton Energy Services, Inc.|Analyzing sand stabilization treatments| CN103206203B|2013-03-28|2015-09-16|重庆大学|The analytical method that the single perforation of oil well shakes out| CN105579664A|2013-08-30|2016-05-11|界标制图有限公司|Reservoir simulator, method and computer program product| WO2016028564A1|2014-08-22|2016-02-25|Schlumberger Canada Limited|Methods for monitoring fluid flow and transport in shale gas reservoirs| GB2548022A|2014-10-24|2017-09-06|Landmark Graphics Corp|Inflow control apparatus, methods,and systems| US10428641B2|2015-04-17|2019-10-01|Landmark Graphics Corporation|Draw-down pressure apparatus, systems, and methods| WO2016191265A1|2015-05-22|2016-12-01|Aramco Services Company|Method for determining unconventional liquid imbibition in low-permeability materials| US10685086B2|2015-09-15|2020-06-16|Conocophillips Company|Avoiding water breakthrough in unconsolidated sands| US10488553B2|2016-04-01|2019-11-26|Baker Hughes, A Ge Company, Llc|Stress tensor computation using Mindlin formulation| CN107894204B|2016-10-04|2020-02-21|财团法人工业技术研究院|Interferometer and imaging method thereof| AU2017382546A1|2016-12-19|2019-07-11|Conocophillips Company|Subsurface modeler workflow and tool| WO2019070252A1|2017-10-04|2019-04-11|Halliburton Energy Services, Inc.|Applying triaxial stresses to a core sample during perforation and flow testing| CN108331574B|2018-01-08|2019-03-05|中国石油大学|A kind of opposite section prediction in debt and the sand control segmentation stage division of shaking out of horizontal well horizontal segment| EP3743747A1|2018-01-24|2020-12-02|Saudi Arabian Oil Company|Hydrocarbon migration and accumulation methods and systems| EP3768939A4|2018-03-21|2021-12-29|Resfrac Corporation|Systems and methods for hydraulic fracture and reservoir simulation| US20210018411A1|2018-08-01|2021-01-21|Halliburton Energy Services, Inc.|Stressed rock perforating-charge testing system| US20210131931A1|2018-12-18|2021-05-06|Halliburton Energy Services, Inc.|Determining when applied stress to a core rock sample has equilibrated in the core rock sample| CN109740210B|2018-12-21|2020-10-13|中国石油大学|Shale gas fracturing underground sand blocking accident real-time risk assessment method and device| CN111425170A|2019-01-08|2020-07-17|中国石油天然气集团有限公司|Multi-objective design method and device for drilling perforation casing based on multi-objective optimization| CN109826622A|2019-03-07|2019-05-31|中国石油大学|The simulation system that simulation sandstone reservoir shakes out| EP3708768A1|2019-03-12|2020-09-16|Chevron U.S.A. Inc.|Coupling a simulator and at least one other simulator| US20210040837A1|2019-08-08|2021-02-11|Saudi Arabian Oil Company|Automated sand grain bridge stability simulator|
法律状态:
2019-01-08| B06F| Objections, documents and/or translations needed after an examination request according art. 34 industrial property law| 2019-07-23| B06U| Preliminary requirement: requests with searches performed by other patent offices: suspension of the patent application procedure| 2020-06-16| B09A| Decision: intention to grant| 2020-09-15| B16A| Patent or certificate of addition of invention granted|Free format text: PRAZO DE VALIDADE: 20 (VINTE) ANOS CONTADOS A PARTIR DE 17/09/2010, OBSERVADAS AS CONDICOES LEGAIS. |
优先权:
[返回顶部]
申请号 | 申请日 | 专利标题 US12/561,886|US8548783B2|2009-09-17|2009-09-17|Computer-implemented systems and methods for controlling sand production in a geomechanical reservoir system| US12/561886|2009-09-17| PCT/US2010/049315|WO2011035146A2|2009-09-17|2010-09-17|Computer-implemented systems and methods for controlling sand production in a geomechanical reservoir system| 相关专利
Sulfonates, polymers, resist compositions and patterning process
Washing machine
Washing machine
Device for fixture finishing and tension adjusting of membrane
Structure for Equipping Band in a Plane Cathode Ray Tube
Process for preparation of 7 alpha-carboxyl 9, 11-epoxy steroids and intermediates useful therein an
国家/地区
|